import re
import sys
import os, os.path
import fcntl
import termios
import struct
import itertools
import copy
import warnings
import datetime
from math import *

# pyfits needed for read/writing ATOMDB FITS files
try:
  import pyfits
except:
  pass
# ==============================================================================
# constants
# ==============================================================================

SPECLETTERS=['S','P','D','F','G','H','I','K','L','M','N','O','Q','R','T','U',\
             'V','W','X','Y','Z']

ATOMSYMBOLS=['A','H','He','Li','Be','B','C','N','O','F','Ne','Na','Mg','Al',\
             'Si','P','S','Cl','Ar','K','Ca','Sc','Ti','V','Cr','Mn','Fe','Co',\
             'Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y','Zr',\
             'Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I',\
             'Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy',\
             'Ho','Er','Tm','Yb','Lu','Hf','Ta','W','Re','Os','Ir','Pt','Au',\
             'Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa','U',\
             'Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No','Lr','Rf','Db',\
             'Sg','Bh','Hs','Mt','Ds','Rg','Uub','Uut','Uuq','Uup','Uuh','Uus',\
             'Uuo']

SUBSHELLS=['0s',\
           '1s',\
           '2s','2p',\
           '3s','3p','3d',\
           '4s','4p','4d','4f',\
           '5s','5p','5d','5f','5g',\
           '6s','6p','6d','6f','6g','6h',\
           '7s','7p','7d','7f','7g','7h','7i']

SSPIECES=[(0,0),\
          (1,0),\
          (2,0),(2,1),\
          (3,0),(3,1),(3,2),\
          (4,0),(4,1),(4,2),(4,3),\
          (5,0),(5,1),(5,2),(5,3),(5,4),\
          (6,0),(6,1),(6,2),(6,3),(6,4),(6,5),\
          (7,0),(7,1),(7,2),(7,3),(7,4),(7,5),(7,6)]

NSHELLS=['A','K','L','M','N','O','P','Q','R','S','T']

ROMAN=['0','I','II','III','IV','V','VI','VII','VIII','IX','X','XI','XII',\
       'XIII','XIV','XV','XVI','XVII','XVIII','XIX','XX','XXI','XXII',\
       'XXIII','XXIV','XXV','XXVI','XXVII','XXVIII','XXIX','XXX','XXXI',\
       'XXXII','XXXIII','XXXIV','XXXV','XXXVI','XXXVII','XXXVIII','XXXIX',\
       'XL']

EQUIV_LS_TERMS={}   # permitted LS terms for equivalent electrons; filled as 
                    # needed by rlf_equiv_terms()

RYDBERG=13.6056923
INVERSECM=109737.3
HC_RYDANG=911.26704
HC_EVNM=1239.8419
GF2A=0.66702

# ==============================================================================
# RAP exception type
# ==============================================================================

class ErrorRAP(Exception):
  def __init__(self,value):
    self.value=value
  def __str__(self):
    return repr(self.value)

# ==============================================================================
# RAP class helper functions
# ==============================================================================

def rlf_filetype(filename):
  """ 
  Check type of file: adf04, adf39l, TCCDW.DAT, olg, NIST, or xstar type 006.
  """
  if '.fits' in filename:
    rlf_have_pyfits(stop=True)
    try:
      a=pyfits.open(filename)
    except:
      pass
    else:
      origin=a[0].header['ORIGIN']
      dtype=a[1].header['HDUCLAS1']
      if origin == 'ATOMDB':
        if dtype == 'E_LEVEL':
          return 'atomdb_ev'
      a.close()

  fin=open(filename)
  lines=[]
  for line in fin:
    lines.append(line)
    if len(lines) == 10: break
  fin.close()

  if '=====' in lines[0]:     # must be TIPTOP
    if 'E(' in lines[1]: return 'tiptop'
    if 'WL(' in lines[1]: return 'tiptop_lines'

  for line in lines:
    if line[:3] == 'FAC':
      itype=int(lines[3].split()[-1])
      if itype == 1:
        return 'fac_lev'
      elif itype == 2:
        return 'fac_tr'
      else:
        break
    if line[:4] == 'SEQ=': return 'adf39l'
    tmp=line.split()
    if len(tmp) == 18 and tmp[-1] == '%' and int(tmp[0]) == 6: return 'x006'
    if len(tmp) in [14,15] and tmp[-1] == '%' and int(tmp[0]) == 50: 
      return 'x050'
    if len(tmp) == 12 and tmp[5] in SPECLETTERS: return 'chianti_elvlc'
    if len(tmp) == 10 and tmp[7] == '-': return 'chianti_wgfa'
    tmp=re.findall(r'(\d+)\s+(\d+)\s+(\d+)\s+(\d+\.)',line)
    if len(tmp) > 0:  return 'adf04'
    tmp=re.findall(r'(LS TERMS)',line)
    if len(tmp) > 0:  return 'tcc'
    if '**********' in line: return 'olg'
    if line.count('|') == 13: return 'nist_lines'
    if line.count('|') in [4,5]: return 'nist'
  raise ErrorRAP('cannot understand file type: '+filename)

# ------------------------------------------------------------------------------

def rlf_loadstates(filename,cup=None,**args):
  """ 
  Read file and create rap_states child object appropriate to detected or
  specified coupling scheme.  If coupling scheme is unspecified then try
  LSJ, LS, then CA.
  """
  if not os.path.exists(filename):
    msg='file, %s, not found'%filename
    raise ErrorRAP(msg)

  # coupling scheme is given...
  if cup == 'CA':
    try:
      out=rap_states_ca(filename,**args)
    except:
      msg='no CA state data found in %s'%filename
    else:
      return out
  if cup == 'LS':
    try:
      out=rap_states_ls(filename,**args)
    except:
      msg='no LS state data found in %s'%filename
    else:
      return out
  if cup == 'LSJ':
    try:
      out=rap_states_lsj(filename,**args)
    except:
      msg='no LSJ state data found in %s'%filename
    else:
      return out
  if cup is not None:
    msg='coupling scheme, %s, not understood'%cup
    raise ErrorRAP(msg)

  # coupling scheme is not given...
  try:
    out=rap_states_lsj(filename,**args)
  except ErrorRAP, e:
    msglsj=e.value
  else:
    return out
  try:
    out=rap_states_ls(filename,**args)
  except ErrorRAP, e:
    msgls=e.value
  else:
    return out
  try:
    out=rap_states_ca(filename,**args)
  except ErrorRAP, e:
    msgca=e.value
  else:
    return out
  msg=[]
  msg.append('LSJ: %s'%msglsj)
  msg.append('LS:  %s'%msgls)
  msg.append('CA:  %s'%msgca)
  msg.append('no state data found of any coupling scheme')
  msg='\n'.join(msg)
  raise ErrorRAP(msg)

# ==============================================================================
# functions acting on ion strings
# ==============================================================================

def rlf_iondecode(ion):
  """ 
  Return tuple, (z,ne), from ion string where z is the nuclear charge and ne is
  the number of electrons.
  """
  # look for sequence notation
  if '-like' in ion:
    try:
      tmp=ion.split('-like ')
      ne=ATOMSYMBOLS.index(tmp[0].strip().capitalize())
      z=ATOMSYMBOLS.index(tmp[1].strip().capitalize())
    except:
      msg="invalid ion string -- sequence notation"
      raise ErrorRAP(msg)
    return (z,ne)

  # look for astrophysical notation
  if ' ' in ion:
    tmp=ion.split()
    try:
      z=ATOMSYMBOLS.index(tmp[0].strip().capitalize())
      q=ROMAN.index(tmp[1].strip())-1
    except:
      msg="invalid ion string -- astrophysical notation"
      raise ErrorRAP(msg)
    else:
      return (z,z-q)
  # look for chemical notation
  tmp=re.findall(r'([A-Za-z][A-Za-z]?)\+?(\d+)\+?',ion)
  if not tmp:
    msg="invalid ion string -- chemical notation"
    raise ErrorRAP(msg)
  z=ATOMSYMBOLS.index(tmp[0][0].capitalize())
  q=int(tmp[0][1])
  return (z,z-q)

# ------------------------------------------------------------------------------

def rlf_ionencode(z,ne=None,form='chemical'):
  """ Take nuclear charge and electron count and return ion string. """
  if ne is None: (z,ne)=z
  q=z-ne
  x=ATOMSYMBOLS[z]
  if form[0].lower() == 'c':
    out=x+str(q)+'+'
  elif form[0].lower() == 'a':
    r=ROMAN[q+1]
    out=x+' '+r
  elif form[0].lower() == 's':
    out="%s-like %s"%(ATOMSYMBOLS[ne],ATOMSYMBOLS[z])
  else:
    msg="invalid output format"
    raise ErrorRAP(msg)
  return out

# ==============================================================================
# functions acting on sub-shells
# ==============================================================================

def rlf_subshell(i):
  """ Return sub-shell string of i'th sub-shell (in standard order). """
  if i < 1: 
    msg='sub-shell index out of range: i=%i'%i
    raise ErrorRAP(msg)
  if i <= len(SUBSHELLS)-1: return SUBSHELLS[i]
  n=int(round(sqrt(2.*float(i)-1.)))
  l=i-1-n*(n-1)/2
  if l > len(SPECLETTERS): 
    msg='angular momenta greater than '+str(len(SPECLETTERS)-1)+\
        ' not supported: i=%i => n,l=%i, %i'%(i,n,l)
    raise ErrorRAP(msg)
  return str(n)+SPECLETTERS[l].lower()

# ------------------------------------------------------------------------------

def rlf_ssindex(ss):
  """ 
  Return the standard order index of subshell, ss. 

  Supports continuum orbitals where the principal quantum number is replaced
  by 'k'.  The returned index is -(l+1) where l is the orbital angular 
  momentum.
  """
  if ss in SUBSHELLS: return SUBSHELLS.index(ss)
  if ss[0].upper() == 'K':                 # check for continuum orbital
    l=SPECLETTERS.index(ss[-1].upper())
    return -(l+1)
  tmp=re.findall(r'(\d+[A-Za-z])',ss)
  n=int(tmp[0][:-1])
  try:
    l=SPECLETTERS.index(tmp[0][-1].upper())
  except:
    msg='illegal symbol for L in spectroscopic notation: %s'%ss
    raise ErrorRAP(msg)
  if l > n-1:
    msg='illegal sub-shell: %s'%ss
    raise ErrorRAP(msg)
  i=n*(n-1)/2+l+1
  return i

# ------------------------------------------------------------------------------

def rlf_ssdecode(ssidx):
  """ Return tuple containing n,l for given sub-shell (string or index). """
  try:
    idx=int(ssidx)
  except:
    idx=rlf_ssindex(ssidx)
  try:
    out=SSPIECES[idx]
  except:
    n=int(round(sqrt(2.*float(idx)-1.)))
    l=idx-1-n*(n-1)/2
    out=(n,l)
  return out

# ------------------------------------------------------------------------------

def rlf_ssencode(n,l,index=False):
  """ Return sub-shell given n and l. """
  ss=str(n)+SPECLETTERS[l].lower()
  if index: ss=rlf_ssindex(ss)
  return ss

# ------------------------------------------------------------------------------

def rlf_ssoccupancy(ssidx):
  """ Return max occupancy of given sub-shell (string or index). """
  (n,l)=rlf_ssdecode(ssidx)
  return 4*l+2

# ------------------------------------------------------------------------------

def rlf_ss2jss(ssidx):
  """ Return list of allowed J subshells given a regular sub-shell. """
  (n,l)=rlf_ssdecode(ssidx)
  if l == 0: return [rlf_jssencode(n,l,+1)]
  return [rlf_jssencode(n,l,-1),rlf_jssencode(n,l,+1)]

# ------------------------------------------------------------------------------

def rlf_sssorter(x,y):
  """ 
  Use this function in the list sort() function in order to arrange items by
  their sub-shell order.
  """
  xx=rlf_ssindex(x)
  yy=rlf_ssindex(y)
  if xx < 0 and yy < 0:
    if xx < yy: return +1
    if xx > yy: return -1
    return 0
  if xx < 0: return +1
  if yy < 0: return -1
  if xx > yy: return +1
  if xx < yy: return -1
  return 0

# ------------------------------------------------------------------------------

def rlf_ss2eissner(idx):
  """ Return symbol for sub-shell in Eissner notation. """
  tidx=idx
  if not isinstance(tidx,int): tidx=rlf_ssindex(tidx)
  if idx <= 0: 
    msg='orbital not supported by Eissner notation; idx: %s'%idx
    raise ErrorRAP(msg)
  if tidx < 10:
    x=str(tidx)
  elif tidx < 36:
    x=chr(tidx-9-1+ord('A'))
  elif tidx < 62:
    x=chr(tidx-35-1+ord('a'))
  else:
    msg='orbital not supported by Eissner notation; idx: %s'%idx
    raise ErrorRAP(msg)
  return x

# ------------------------------------------------------------------------------

def rlf_eissner2ss(x,index=False):
  """ Convert sub-shell symbol from Eissner notation into standard form. """
  x=str(x)
  if len(x) > 1: 
    msg='Eissner orbital symbol too long: %s'%x
    raise ErrorRAP(msg)
  if x >= '1' and x <= '9':
    ss=int(x)
  elif x >= 'A' and x <= 'Z':
    ss=ord(x)-ord('A')+1+9
  elif x >= 'a' and x <= 'z':
    ss=ord(x)-ord('a')+1+35
  else:
    msg='invalid symbol for orbital in Eissner notation: %s'%x
    raise ErrorRAP(msg)
  if index: return ss
  return rlf_subshell(ss)

# ==============================================================================
# functions acting on JJ-coupled sub-shells
# ==============================================================================

def rlf_jsubshell(i):
  n=int(sqrt(float(i)-0.5))+1
  d=i-(n-1)**2
  l=d/2
  s=d%2
  out="%i%s"%(n,SPECLETTERS[l].lower())
  if s == 0:
    out+='-'
  else:
    out+='+'
  return out

# ------------------------------------------------------------------------------

def rlf_jssindex(ss):
  (n,l,s)=rlf_jssdecode(ss)
  return (n-1)**2+2*l+(s+1)/2

# ------------------------------------------------------------------------------

def rlf_jssdecode(iorss):
  if isinstance(iorss,int): iorss=rlf_jsubshell(iorss)
  try:
    n=int(iorss[:-2])
    l=SPECLETTERS.index(iorss[-2].upper())
    if '+' in iorss: 
      s=1
    elif '-' in iorss:
      s=-1
    else:
      raise
  except:
    msg='cannot understand J-subshell'
    raise ErrorRAP(msg)
  return (n,l,s)

# ------------------------------------------------------------------------------

def rlf_jssencode(n,l,s,index=False):
  if s not in [-1,+1]:
    msg='invalid value for sigma for J-subshell'
    raise ErrorRAP(msg)
  if l == 0 and s == -1:
    msg='sigma must be +1 when l=0'
    raise ErrorRAP(msg)
  idx=(n-1)**2+2*l+(s+1)/2
  if index: return idx
  return rlf_jsubshell(idx)

# ------------------------------------------------------------------------------

def rlf_jssoccupancy(iorss):
  (n,l,s)=rlf_jssdecode(iorss)
  return 2*l+s+1

# ------------------------------------------------------------------------------

def rlf_jsssorter(x,y):
  xx=rlf_jssindex(x)
  yy=rlf_jssindex(y)
  return cmp(xx,yy)

# ==============================================================================
# functions acting term strings
# ==============================================================================

def rlf_termsorter(x,y):
  """
  Use this function in the list sort() in order to arrange entries according 
  to term order, e.g. tlist.sort(rlf_termsorter).
  """
  return cmp(rlf_termdecode(x),rlf_termdecode(y))

# ------------------------------------------------------------------------------

def rlf_termdecode(term):
  """ Convert term string to tuple of 2S+1 and L. """
  try:
    g=int(term[:-1])
    l=SPECLETTERS.index(term[-1])
  except:
    msg='bad term string: %s'%term
    raise ErrorRAP(msg)
  return (g,l)

# ------------------------------------------------------------------------------

def rlf_termencode(g,l):
  """ Return term string given 2S+1 and L. """
  if l < 0 or l > len(SPECLETTERS)-1:
    msg='angular momentum out of supported range; L=%s'%l
    raise ErrorRAP(msg)
  return str(g)+SPECLETTERS[l]

# ------------------------------------------------------------------------------

def rlf_termweight(term):
  """ Return term weight: (2S+1)(2L+1). """
  (g,l)=rlf_termdecode(term)
  return g*(2*l+1)

# ------------------------------------------------------------------------------

def rlf_term_2jvals(term):
  """ Return list of allowed 2J values for given term. """
  (g,l)=rlf_termdecode(term)
  j2min=abs(g-(2*l+1))
  j2max=g+(2*l+1)-2
  return range(j2min,j2max+1,2)

# ------------------------------------------------------------------------------

def _rlf_calc_equiv_lsterms(lorb,ne):
  """ 
  Calculate dictionary of LS terms of orbital with angular momentum, lorb, 
  and having ne electrons.

  The returned dictionary has keys composed of tuples having the form 
  (2S+1,L,p), where 2S+1,L is the usual term identifiers and p is the parity.
  The value of the dictionary with a given key is the number of terms in that
  list with that key.
  """
  def to_binary(x,npad=0):
    """ Convert integer into binary representation via list """
    xx=int(x)
    out=[]
    while xx > 0:
      (xx,y)=divmod(xx,2)
      out.append(y)
    if npad > len(out):
      for i in range(len(out),npad): out.append(0)
    out.reverse()
    return out
  
  def array_mult(x,y):
    """ Multiply elements of two arrays together to form new array """
    if len(x) != len(y): raise ErrorRAP('arrays must be of equal size')
    out=[]
    for i in range(len(x)):
      out.append(x[i]*y[i])
    return out
  
  def dot_prod(x,y):
    """ Take dot-product of two arrays of equal size """
    vec=array_mult(x,y)
    return sum(vec)

  # short-cut for single electron (hole)
  if ne == 1 or ne == 4*lorb+1:
    termlist={}
    key=(2,lorb,(-1)**lorb)
    termlist[key]=1
    return termlist

  # define orbital
  ll=lorb
  nelec=ne
  parity=(1-(-1)**(ll*nelec))/2
  
  # form arrays defining angular momenta and spins values
  lval=[]
  for i in range(2*ll+1):
    lval.append(-ll+i)
  lval*=2
  sval=[-1]*(2*ll+1)+[+1]*(2*ll+1)
  
  # calculate all possible combinations of L and S
  nbox=4*ll+2
  nmax=2**nbox-1
  stot=[]
  ltot=[]
  for n in range(1,nmax+1):
    arr=to_binary(n,npad=nbox)
    if sum(arr) == nelec:
      stot.append(dot_prod(sval,arr))
      ltot.append(dot_prod(lval,arr))

  # form matrix with frequency distribution of L and S
  lmax=max(ltot)
  smax=max(stot)
  col=[]
  l=-lmax
  while l <= lmax:
    col.append(0)
    l+=1
  table=[]
  s=-smax
  while s <= smax:
    table.append(col[:])
    s+=2
  for i in range(len(ltot)):
    jl=ltot[i]+lmax
    js=(stot[i]+smax)/2
    table[js][jl]+=1

  # get list of terms
  termlist={}
  jlmax=(len(col)+1)/2-1
  jsmax=(len(table)+1)/2-1
  js=0
  while js <= jsmax:
    jl=0
    s=2*js-smax
    while jl <= jlmax:
      l=jl-lmax
      n=table[js][jl]
      if n > 0:
        key=(abs(s)+1,abs(l),parity)
        for kl in range(jl,len(col)-jl):
          for ks in range(js,len(table)-js):
            table[ks][kl]=table[ks][kl]-n
        termlist[key]=n
      jl=jl+1
    js=js+1
  return termlist

# ------------------------------------------------------------------------------

def _rlf_lsterm_couple(t1,t2,numterms=1):
  """ 
  Internal function for coupling two terms together; each term is a tuple of
  the form, (2S+1,L,p).  Returns a LS term list with the number of terms set
  to numterms.
  """
  (s1,l1,p1)=t1
  (s2,l2,p2)=t2
  n=int(numterms)
  smin=abs(s1-s2)
  smax=abs(s1+s2-2)
  lmin=abs(l1-l2)
  lmax=abs(l1+l2)
  p=p1+p2
  if p == 2: p=0
  termlist={}
  for s in range(smin,smax+1,2):
    for l in range(lmin,lmax+1):
      nt=(s+1,l,p)
      termlist[nt]=n
  return termlist

# ------------------------------------------------------------------------------

def rlf_equiv_lsterms(lorb,ne):
  """ 
  Return dictionary of LS terms from an orbital with angular momentum, lorb, 
  having ne equivalent electrons.

  The dictionary of terms is only calculated once and then stored in the 
  EQUIV_LS_TERMS global variable (reset after session ends).
  """
  try:
    ll=int(lorb)
    nn=int(ne)
  except:
    return False
  if nn > 4*ll+2: 
    msg='too many electrons for shell: lorb,ne: %s, %s'%(lorb,ne)
    raise ErrorRAP(msg)
  ss=(ll,nn)
  if ss in EQUIV_LS_TERMS.keys(): return EQUIV_LS_TERMS[ss]
  termlist=_rlf_calc_equiv_lsterms(ll,nn)
  EQUIV_LS_TERMS[ss]=termlist
  return termlist

# ------------------------------------------------------------------------------

def rlf_lsterm_merge(termlist1,termlist2):
  """ Merge termlists together and return result. """
  if not isinstance(termlist1,dict):
    msg='input term lists must be dictionaries, received: %s'%type(termlist1)
    raise ErrorRAP(msg)
  if not isinstance(termlist2,dict):
    msg='input term lists must be dictionaries, received: %s'%type(termlist2)
    raise ErrorRAP(msg)
  termlist=termlist1.copy()
  for t,n in termlist2.iteritems():
    if t not in termlist:
      termlist[t]=n
    else:
      termlist[t]+=n
  return termlist

# ------------------------------------------------------------------------------

def rlf_lsterm_couple(termlist1,termlist2):
  """ 
  Couple all terms from termlist1 with terms from termlist2 and return the
  resulting dictionary.
  """
  if not isinstance(termlist1,dict):
    msg='input term lists must be dictionaries, received: %s'%type(termlist1)
    raise ErrorRAP(msg)
  if not isinstance(termlist2,dict):
    msg='input term lists must be dictionaries, received: %s'%type(termlist2)
    raise ErrorRAP(msg)
  termlist={}
  for t1,n1 in termlist1.iteritems():
    for t2,n2 in termlist2.iteritems():
      n=n1*n2
      termlist=rlf_lsterm_merge(termlist,_rlf_lsterm_couple(t1,t2,n))
  return termlist

# ------------------------------------------------------------------------------

def rlf_num_lsterms(termlist):
  """ Return number of terms in list """
  if not isinstance(termlist,dict):
    msg='input term list must be a dictionary, received: %s'%type(termlist)
    raise ErrorRAP(msg)
  total=0
  for n in termlist.itervalues(): 
    total+=n
  return total

# ==============================================================================
# general configuration functions
# ==============================================================================

def rlf_confsorter(x,y):
  """
  Use this function in the list sort() in order to arrange entries according 
  to configuration order, e.g. alist.sort(rlf_confsorter).
  """
  if not isinstance(x,rap_config):
    msg='only rap_config types can be sorted; received: %s'%type(x)
    raise ErrorRAP(msg)
  if not isinstance(y,rap_config):
    msg='only rap_config types can be sorted; received: %s'%type(y)
    raise ErrorRAP(msg)
  ssx=x.subshells(indices=True)
  ssy=y.subshells(indices=True)
  if x.num_electrons() < y.num_electrons(): return -1
  if x.num_electrons() > y.num_electrons(): return +1
  ssmax=max(max(ssx),max(ssy))
  for i in range(1,ssmax+1):
    if i not in ssx and i not in ssy: continue
    if i not in ssx: return +1
    if i not in ssy: return -1
    wx=x.occupancy(i)
    wy=y.occupancy(i)
    if wx < wy: return +1
    if wx > wy: return -1
    if i == ssmax: return 0
  raise ErrorRAP('sanity check - should have returned a value')

# ------------------------------------------------------------------------------

def rlf_jconfsorter(x,y):
  """
  Use this function in the list sort() in order to arrange entries according 
  to configuration order, e.g. alist.sort(rlf_confsorter).
  """
  if not isinstance(x,rap_jconfig):
    msg='only rap_jconfig types can be sorted; received: %s'%type(x)
    raise ErrorRAP(msg)
  if not isinstance(y,rap_jconfig):
    msg='only rap_jconfig types can be sorted; received: %s'%type(y)
    raise ErrorRAP(msg)
  ssx=x.subshells(indices=True)
  ssy=y.subshells(indices=True)
  if x.num_electrons() < y.num_electrons(): return -1
  if x.num_electrons() > y.num_electrons(): return +1
  ssmax=max(max(ssx),max(ssy))
  for i in range(1,ssmax+1):
    if i not in ssx and i not in ssy: continue
    if i not in ssx: return +1
    if i not in ssy: return -1
    wx=x.occupancy(i)
    wy=y.occupancy(i)
    if wx < wy: return +1
    if wx > wy: return -1
    if i == ssmax: return 0
  raise ErrorRAP('sanity check - should have returned a value')

# ------------------------------------------------------------------------------

def rlf_genconf(aconf,force=False):
  """
  Generate list of configurations from a configuration string which includes
  n- and l- variables.  

  Example: '2l5 nl n=3,4 l=0-2' will expand into a 
  list of 18 configurations.  Multiple variables can be used as long as the
  n-variables begin with 'n' and the l-variables begin with 'l'.  Also, the
  variables cannot contain any numbers.  Finally, to avoid confusion, the 
  n-variable cannot contain any l's.  For each n-variable, a range must be
  specified (as in the above example).  The l-variables will take all values
  by default unless a range is given.  For cases where an l-electron is
  permissible (L=8 electrons), the routine will assume you want the single
  electron unless the force variable is True or 'l' is used as a variable 
  somewhere else in the string.  Note that no spaces are allowed in the 
  range definitions (c.f. n = 3,4 is legal; n=3, 4 is illegal).  

  Examples:  
    2l5                          ( 3 configurations)
    3s 3l                        ( 3 configurations)
    nl npl n=3,4 np=3-5 l=0-2    (54 configurations)
  """
  def fill_conf(n,lvals,ne,conf0=''):
  # A recursive function to generate a list of configurations with ne 
  # electrons for a given n-value for values of l in lvals.  The variable
  # conf0 is used in the recursive calling process.
    ll=lvals[0]
    if len(lvals) == 1:
      if ne > 4*ll+2: return None
      if SPECLETTERS[ll].lower() == 'l':
        ss=str(n)+'LTEMP'+str(ne)
      else:
        ss=str(n)+SPECLETTERS[ll].lower()+str(ne)
      return [(conf0+' '+ss).strip()]
    else:
      out=[]
      for k in range(ne+1):
        tne=ne-k
        if tne > 4*ll+2: continue
        if SPECLETTERS[ll].lower() == 'l':
          ss=str(n)+'LTEMP'+str(tne)
        else:
          ss=str(n)+SPECLETTERS[ll].lower()+str(tne)
        tconf=(conf0+' '+ss).strip()
        tmp=fill_conf(n,lvals[1:],k,tconf)
        if isinstance(tmp,list):
          out.extend(tmp)
        elif tmp:
          out.append(tmp.strip())
      return out

  def parse_range(x):
  # Return a list of values from the user input range of n or l.
  # (Example: 3,5-7 -> [3,5,6,7] )
    tmp1=re.findall(r'^\d+',x)
    tmp2=re.findall(r'([,-])(\d+)',x)
    if len(tmp1) == 0: return False
    s=[tmp1[0]]
    gotdash=False
    if len(tmp2) != 0:
      for t in tmp2:
        if t[0] == ',':
          s.append(t[1])
          gotdash=False
        else:
          if gotdash == True: return False
          s[-1]=str(s[-1])+t[0]+t[1]
          gotdash=True
    out=[]
    for it in s:
      if '-' not in it:
        out.append(int(it))
      else:
        tmp=re.findall(r'(\d+)',it)
        out.extend(range(int(tmp[0]),int(tmp[1])+1))
    return out
  
  # count number of subshells
  tmp=re.findall(r'([^\s]+)',str(aconf))
  sscount=0
  for i in range(len(tmp)):
    t=tmp[i]
    if '=' not in t: 
      sscount=sscount+1
      try:
        ne=int(t[-1])
      except:
        tmp[i]=t+'1'
  conf=''
  for t in tmp: conf=conf+' '+t
  conf=conf.strip()

  # identify metavariables
  tnmeta=[]
  tlmeta=[]
  ce=re.findall(r'(n[A-Za-km-z]*|\d+)(l[A-Za-z]*|[a-z])(\d+)',conf)
  if len(ce) == 0: 
    msg='invalid configuration: %s'%conf
    raise ErrorRAP(msg)
  for ss in ce:
    # check if 'n' is an integer
    try:
      n=int(ss[0])
    except:
      if ss[0] not in tnmeta: tnmeta.append(ss[0])
      n=0
    # check if 'l' is spectroscopic letter
    if ss[1][0] != 'l':
      if ss[1].upper() not in SPECLETTERS: 
        msg="l metavariable must start with 'l': %s"%ss
        raise ErrorRAP(msg)
      l=SPECLETTERS.index(ss[1].upper())
    else:
      l=-1
      if len(ss[1]) > 1:
        if ss[1] not in tlmeta: tlmeta.append(ss[1])
      # if l='l', check if it is legal for given n
      # always metavar if force is True
      else:
        if n < 9 or force == True: 
          if ss[1] not in tlmeta: tlmeta.append(ss[1])
    # check if legal number of electrons
    ne=int(ss[2])
    if l >= 0:
      if ne > 4*l+2: 
        msg='too many electrons for l=%i'%l
        raise ErrorRAP(msg)
    elif n > 0:
      if ne > 2*n*n: 
        msg='too many electrons for n=%i'%n
        raise ErrorRAP(msg)
  # check that all subshells were legal
  if len(ce) != sscount:
    msg='invalid configuration string: %s'%conf
    raise ErrorRAP(msg)

  # parse range information
  nmeta={}
  lmeta={}
  tmp=re.findall(r'([^\s]+)=([^\s]+)',conf)
  for t in tmp:
    if t[0] in tnmeta:
      r=parse_range(t[1])
      if r != False: nmeta[t[0]]=r
    elif t[0] in tlmeta:
      r=parse_range(t[1])
      if r != False: lmeta[t[0]]=r
    else:
      continue
  for n in tnmeta:
    if n not in nmeta:
      msg='must specify range for all n metavariables'
      raise ErrorRAP(msg)
  for l in tlmeta:
    if l not in lmeta:
      lmeta[l]=[]

  # get list of possibilities for each subshell
  nmin=None
  lmin=None
  poss=[]
  null=[]
  for ss in ce:
    if ss[0] in nmeta:
      nval=nmeta[ss[0]]
    else:
      nval=[int(ss[0])]
    if nmin is None: nmin=min(nval)
    if ss[1] in lmeta:
      lval=lmeta[ss[1]]
    else:
      lval=[SPECLETTERS.index(ss[1].upper())]
    tlmin=min(lval,0)
    if lmin is None: lmin=tlmin
    if min(nval) == nmin:
      if tlmin < lmin: lmin=tlmin
    elif min(nval) < nmin:
      lmin=tlmin
      nmin=min(nval)
    ne=ss[2]
    p=[]
    for n in nval:
      if len(lval) == 0:
        lrange=range(int(n))
      else:
        lrange=[]
        for la in lval:
          if la < int(n): lrange.append(la)
      if(int(ne)) == 0:
        for l in lrange: null.append(str(n)+SPECLETTERS[l].lower())
      p.extend(fill_conf(int(n),lrange,int(ne)))
    poss.append(p)

  # remove temporary variable LTEMP
  for i in range(len(poss)):
    for j in range(len(poss[i])):
      try:
        poss[i][j]=poss[i][j].replace('LTEMP','l')
      except:
        pass

  # construct configurations
  out=poss.pop(0)
  for p in poss:
    tmp=[]
    for ap in p:
      for ss in out:
        tmp.append(ss+' '+ap)
    out=tmp

  # form and sort configuration
  nout=[]
  ne=None
  for i in range(len(out)): 
    try:
      c=rap_config(out[i],ne=ne)
    except ErrorRAP, e:
      if 'maximum occupancy exceeded' in e.value: continue
    nout.append(c)
    if ne is None: ne=nout[-1].num_electrons()
  nout.sort(rlf_confsorter)
  return nout

# ==============================================================================
# angular algebra
# ==============================================================================

def pyang_triangle(j1,j2,j3):
  """ Return True or False if triangle relations are met. """
  if j1+j2 < j3: return False
  if j2+j3 < j1: return False
  if j3+j1 < j2: return False
  return True

# ------------------------------------------------------------------------------

def pyang_wigner3j(j1,j2,j3,m1,m2,m3):
  """ Return Wigner 3-j symbol. """
  if m1+m2+m3 != 0.0: return 0.0
  if int(2*(j1+j2+j3)) % 2 != 0: return 0.0
  if not pyang_triangle(j1,j2,j3): return 0.0

  out=factorial(j1+j2-j3)*factorial(j1-j2+j3)*factorial(-j1+j2+j3)
  out*=factorial(j1-m1)*factorial(j1+m1)
  out*=factorial(j2-m2)*factorial(j2+m2)
  out*=factorial(j3-m3)*factorial(j3+m3)
  out=float(out)/float(factorial(j1+j2+j3+1))
  out=sqrt(out)

  kmin=int(max(0,j2-j3-m1,j1-j3+m2))
  kmax=int(min(j1+j2-j3,j1-m1,j2+m2))
  ksum=0.0
  for k in range(kmin,kmax+1):
    t=factorial(k)*factorial(j1+j2-j3-k)*factorial(j1-m1-k)
    t*=factorial(j2+m2-k)*factorial(j3-j2+m1+k)*factorial(j3-j1-m2+k)
    ksum+=(-1.)**k/float(t)

  out*=(-1.)**int(j1-j2-m3)*ksum
  return out

# ------------------------------------------------------------------------------

def pyang_wigner6j(j1,j2,j3,l1,l2,l3):
  """ Return Wigner 6-j symbol. """

  def tri(a,b,c):
    if int(2*(a+b+c)) % 2 != 0: return 0.0
    if not pyang_triangle(a,b,c): return 0.0
    out=factorial(a+b-c)*factorial(a-b+c)*factorial(-a+b+c)
    out=float(out)/float(factorial(a+b+c+1.))
    out=sqrt(out)
    return out

  out=tri(j1,j2,j3)*tri(j1,l2,l3)*tri(l1,j2,l3)*tri(l1,l2,j3)
  if out == 0.0: return out

  kmin=int(max(j1+j2+j3,j1+l2+l3,l1+j2+l3,l1+l2+j3))
  kmax=int(min(j1+j2+l1+l2,j2+j3+l2+l3,j3+j1+l3+l1))
  ksum=0.0
  for k in range(kmin,kmax+1):
    t=factorial(k-j1-j2-j3)*factorial(k-j1-l2-l3)*factorial(k-l1-j2-l3)
    t*=factorial(k-l1-l2-j3)*factorial(j1+j2+l1+l2-k)
    t*=factorial(j2+j3+l2+l3-k)*factorial(j3+j1+l3+l1-k)
    ksum+=(-1.)**k*float(factorial(k+1))/float(t)

  out*=ksum
  return out

# ------------------------------------------------------------------------------

def pyang_wigner9j(j1,j2,j3,j4,j5,j6,j7,j8,j9):
  """ Return Wigner 9-j symbol. """
  j2min=int(2.*max(abs(j1-j9),abs(j4-j8),abs(j2-j6)))
  j2max=int(2.*max(j1+j9,j4+j8,j2+j6))
  jsum=0.0
  for jj2 in range(j2min,j2max+1,2):
    j=0.5*float(jj2)
    t=(-1.)**(jj2)*float(jj2+1)*pyang_wigner6j(j1,j4,j7,j8,j9,j)
    t*=pyang_wigner6j(j2,j5,j8,j4,j,j6)*pyang_wigner6j(j3,j6,j9,j,j1,j2)
    jsum+=t
  return jsum

# ==============================================================================
# numerical functions
# ==============================================================================

def rlf_aitken(x,y,x0,logx=False,logy=False):
    """
    Use aitken method to interpolate data points, (x_i,y_i) in input arrays x 
    and y at point x0.
    """
    if not isinstance(x,list) and not isinstance(x,tuple):
      msg='input arrays must be lists; received: %s'%type(x)
      raise ErrorRAP(msg)
    if not isinstance(y,list) and not isinstance(y,tuple):
      msg='input arrays must be lists; received: %s'%type(y)
      raise ErrorRAP(msg)
    if len(x) != len(y):
      msg='input arrays have different lengths: len(x),len(y)=%s,%s'% \
          (len(x),len(y))
      raise ErrorRAP(msg)
    if logx:
        tx=[]
        for xx in x: tx.append(log10(xx))
        tx0=log10(x0)
    else:
        tx=x[:]
        tx0=x0
    if logy:
        ty=[]
        for yy in y: ty.append(log10(yy))
    else:
        ty=y[:]
    npts=len(tx)
    h=[]
    yt=[]
    for k in range(npts):
        h.append(tx0-tx[k])
        yt.append(ty[k])
    for k in range(npts-1):
        for l in range(k+1,npts):
            yt[l]=(h[k]*yt[l]-h[l]*yt[k])/(tx[l]-tx[k])
    out=yt[npts-1]
    if logy: out=10.0**out
    return out

# ==============================================================================
# xstar functions
# ==============================================================================

def xstar_ionencode(zeff,ne):
  """ Return xstar ion number. """
  if zeff < 1:
    msg='nuclear charge not physical: Z=%s'%zeff
    raise ErrorRAP(msg)
  if ne > zeff:
    msg='negative ions not supported: Z,ne:%s,%s'%(zeff,ne)
    raise ErrorRAP(msg)
  return zeff*(zeff+1)/2+1-ne

# ------------------------------------------------------------------------------

def xstar_iondecode(ionno):
  """ Return tuple, (Z,ne), for the given xstar ion number. """
  if ionno < 1:
    msg='xstar ion number must be positive: ionno=%s'%ionno
    raise ErrorRAP(msg)
  zeff=int(ceil(-0.5+0.5*sqrt(8.*float(ionno)-3.)))
  ne=zeff*(zeff+1)/2+1-ionno
  return (zeff,ne)

# ==============================================================================
# ADAS functions
# ==============================================================================

def rlf_adas2float(x):
  """ 
  Takes a number in an ADAS file (no 'E' given) and gives a floating point 
  number.
  """
  return float(x[:-3]+'E'+x[-3:])

# ==============================================================================
# functions acting AUTOSTRUCTURE olg files
# ==============================================================================

def rlf_olg_params(filename,indices=False,terse=False):
  """
  Return a dictionary of Thomas-Fermi scaling parameters in the given 
  AUTOSTRUCTURE output file where the key is the subshell indentifier, eg. 2p.
  If terse then return a list of scaling parameters with no indentifier.
  """
  if rlf_filetype(filename) != 'olg': 
    msg="invalid file type for %s"%filename
    raise ErrorRAP(msg)
  fin=open(filename,'r')
  while 1:
    line=fin.readline()
    if line == '': 
      msg="end-of-file reached before orbital specification"
      raise ErrorRAP(msg)
    if 'INTERNAL ORBITALS' in line: break
  torbs=[]
  while 1:
    line=fin.readline()
    if 'CONFIGURATION' in line: break
    if 'N L' in line: torbs.extend(line.split()[2:])
  norbs=len(torbs)/2
  orbs=[]
  for i in range(norbs):
    ss=str(torbs[2*i])+SPECLETTERS[int(torbs[2*i+1])].lower()
    if indices: ss=rlf_ssindex(ss)
    orbs.append(ss)
  while 1:
    line=fin.readline()
    if line == '':
      msg="end-of-file reached before scaling parameters"
      raise ErrorRAP(msg)
    if 'SCALING' in line: break
  bline=line[:-1]
  while 1:
    line=fin.readline()
    if line == '' or 'I' in line: break
    bline+=line[:-1]
  tmp=re.findall(r'(\d+\.\d+)',bline)
  params=[]
  for t in tmp: params.append(float(t))
  fin.close()
  if terse: return params
  if len(orbs) != len(params):
    msg="number of scaling parameters does not match number of orbitals"
    raise ErrorRAP(msg)
  return dict(zip(orbs,params))

# ==============================================================================
# compatibility functions
# ==============================================================================

def rlf_have_pyfits(stop=False):
  """ Return True if the pyfits module is installed. """
  if 'pyfits' in sys.modules.keys(): return True
  if stop:
    msg='module, pyfits, is required'
    raise ErrorRAP(msg)
  return False

# ------------------------------------------------------------------------------

# implement sorted() function for Python v2.3
try:
  sorted
except NameError: # doesn't have "sorted"
  def sorted(arr):
    out=arr[:]
    out.sort()
    return out

# ------------------------------------------------------------------------------

# implement reversed() function for Python v2.3
try:
  reversed
except NameError: # doesn't have "reversed"
  def reversed(arr):
    out=arr[:]
    out.reverse()
    return out

# ------------------------------------------------------------------------------

# implement itertools.product
try:
  iter_product=itertools.product
except:
  def iter_product(*args, **kwds):
    # product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
    # product(range(2), repeat=3) --> 000 001 010 011 100 101 110 111
    pools = map(tuple, args) * kwds.get('repeat', 1)
    result = [[]]
    for pool in pools:
      result = [x+[y] for x in result for y in pool]
    for prod in result:
      yield tuple(prod)

# ------------------------------------------------------------------------------

# implement itertools.combinations
try:
  iter_combinations=itertools.combinations
except:
  def iter_combinations(iterable, r):
    # combinations('ABCD', 2) --> AB AC AD BC BD CD
    # combinations(range(4), 3) --> 012 013 023 123
    pool = tuple(iterable)
    n = len(pool)
    if r > n:
      return
    indices = range(r)
    yield tuple([pool[i] for i in indices])
    while True:
      for i in reversed(range(r)):
        if indices[i] != i + n - r:
          break
      else:
        return
      indices[i] += 1
      for j in range(i+1, r):
        indices[j] = indices[j-1] + 1
      yield tuple([pool[i] for i in indices])

# ------------------------------------------------------------------------------

# implement itertools.combinations_with_replacement
try:
  iter_combinations_with_replacement=itertools.combinations_with_replacement
except:
  def iter_combinations_with_replacement(iterable, r):
    """ Python 2.7 function in itertools. """
    pool = tuple(iterable)
    n = len(pool)
    for indices in iter_product(range(n), repeat=r):
      if sorted(indices) == list(indices):
        yield tuple([pool[i] for i in indices])



# ==============================================================================
# miscellaneous functions
# ==============================================================================

def rlf_termsize():
  """
  Return tuple containing the height and width of the terminal window.
  """
  return struct.unpack('hh',fcntl.ioctl(0, termios.TIOCGWINSZ, '1234'))

# ------------------------------------------------------------------------------

def rlf_equiv(alist):
  """ Return True if all items in list, alist, are identical. """
  if alist.count(alist[0]) == len(alist): return True
  return False

# ------------------------------------------------------------------------------

def rlf_binom(a,b):
  """ Return binomial of (a b). """
  if b > a: raise ErrorRAP('invalid binomial coefficients: b > a')
  tb=max(b,a-b)
  top=1.0
  for i in range(tb+1,a+1): top*=i
  bot=1.0
  for i in range(1,a-tb+1): bot*=i
  return top/bot

# ------------------------------------------------------------------------------

def rlf_ordinal(i):
  """ Return ordinal of given integer, e.g. 1 yields 1st. """
  if i%10 == 1 and i%100 != 11:
    return str(i)+'st'
  elif i%10 == 2 and i%100 != 12:
    return str(i)+'nd'
  elif i%10 == 3 and i%100 != 13:
    return str(i)+'rd'
  else:
    return str(i)+'th'

# ------------------------------------------------------------------------------

def rlf_wraptext(text,length,indent=0):
  """
  Take string and return a list of strings where the length of each is at most
  the given length.  Separation of the initial string will only occur at 
  spaces.  An optional indentation can be given for all lines following the
  first.
  """
  if len(text) < length: return text
  out=[]
  pieces=text.split()
  line=pieces.pop(0)
  while len(pieces) > 0:
    if len(line)+len(pieces[0]) > length:
      out.append(line)
      line=" "*indent+pieces.pop(0)
    else:
      line+=" "+pieces.pop(0)
  out.append(line)
  return out

################################################################################
class rap_states:
  """ Superclass for storing state data: CA, LS, or IC. """
  def __init__(self):
    self.clear()

# ------------------------------------------------------------------------------

  def __len__(self):
    return self.num_levels()

# ------------------------------------------------------------------------------

  def __getitem__(self,key):
    return self.label(key)

# ------------------------------------------------------------------------------

  def __contains__(self,i):
    return self.v_energy.__contains__(i)

# ------------------------------------------------------------------------------

  def __iter__(self):
    return self.v_energy.__iter__()

# ------------------------------------------------------------------------------

  def _mapme(self,idx,imap=True):
    tidx=idx
    if imap: tidx=self.lmap(tidx)
    if tidx is None: return None
    if tidx not in self: raise ErrorRAP('level index out of range')
    return tidx

# ------------------------------------------------------------------------------

  def _mapme_r(self,idx,imap=True):
    tidx=idx
    if imap: tidx=self.lmap_reverse(tidx)
    if tidx is None: return None
    return tidx

# ------------------------------------------------------------------------------

  def _map_straight(self,other,fromlist,tolist):
    """ 
    General mapping function based on energy order for levels in fromlist and
    tolist grouped by symmetry.
    """
    if type(self) != type(other):
      msg="cannot map different types"
      raise ErrorRAP(msg)

    out={}
    for key,ifrom in fromlist.iteritems():
      if key not in tolist: continue
      ito=tolist[key]

      # sort indices by energy
      tmp=sorted([(self.energy(i),i) for i in ifrom])
      tifrom=[j for (i,j) in tmp]
      tmp=sorted([(other.energy(i),i) for i in ito])
      tito=[j for (i,j) in tmp]

      tmap={}
      while len(tifrom) > 0 and len(tito) > 0: 
        tmap[tito.pop(0)]=tifrom.pop(0)
      for i in ito: 
        if i not in tmap: tmap[i]=None
      for k,v in tmap.iteritems(): out[k]=v
    return out

# ------------------------------------------------------------------------------

  def _map_by_energy(self,other,fromlist,tolist):
    """ 
    General mapping function minimizing energy differences for levels in 
    fromlist and tolist grouped by symmetry.
    """
    if type(self) != type(other):
      msg="cannot map different types"
      raise ErrorRAP(msg)

    out={}
    for key,ifrom in fromlist.iteritems():
      if key not in tolist: continue
      ito=tolist[key]

      # sort indices by energy
      tmp=sorted([(self.energy(i),i) for i in ifrom])
      tifrom=[j for (i,j) in tmp]
      tmp=sorted([(other.energy(i),i) for i in ito])
      tito=[j for (i,j) in tmp]

      tmap={}
      for i in tifrom:
        if len(tito) == 0: continue
        ediff=[abs(self.energy(i)-other.energy(j)) for j in tito]
        jj=ediff.index(min(ediff))
        tmap[tito[jj]]=i
        tito=tito[jj:]
      for i in ito: 
        if i not in tmap: tmap[i]=None
      for k,v in tmap.iteritems(): out[k]=v
    return out

# ------------------------------------------------------------------------------

  def clear(self):
    """ Erase all data. """
    self.v_cup=None            # coupling scheme defined by child type
    self.v_zeff=None
    self.v_nelec=None
    self.v_core_energy=None    # in Ry
    self.v_ion_thresh=None     # in Ry
    self.v_parity={}
    self.v_energy={}           # in Ry
    self.v_lmap={}
    self.v_lmap_r={}           # reverse map

# ------------------------------------------------------------------------------

  def clear_levels(self):
    """ Just erase level data. """
    self.v_parity={}
    self.v_energy={}
    self.v_lmap={}
    self.v_lmap_r={}

# ------------------------------------------------------------------------------

  def copy(self):
    """ Return a copy. """
    return copy.deepcopy(self)

# ------------------------------------------------------------------------------

  def loadfile(self,filename):
    """ Load in level information from level data file. """
    raise 'No input data types defined'

# ------------------------------------------------------------------------------

  def coupling(self):
    """ Return coupling scheme. """
    return self.v_cup

# ------------------------------------------------------------------------------

  def indices(self,imap=False):
    """ Return list of level indices. """
    out=self.v_energy.keys()
    if imap is True: out=map(self._mapme_r,out)
    out.sort()
    return out

# ------------------------------------------------------------------------------

  def reindex(self,remap=None):
    """ 
    Reassign level indices.  Indices are remapped according to dictionary,
    remap, old (key) -> new (value).  If remap is None, then reindex according
    to energy order.
    """
    if remap:
      err=False
      for i in self.indices():
        if i not in remap:
          err=True
          break
      if err:
        raise ErrorRAP("all indices must be represented using reindex")
      imap=remap
    else:
      energies=[]
      for i,e in self.v_energy.iteritems(): energies.append((e,i))
      energies.sort()
      imap={}
      for j in range(len(energies)): imap[energies[j][1]]=j+1
    self.v_parity=dict([(imap[i],val) for (i,val) in self.v_parity.items()])
    self.v_energy=dict([(imap[i],val) for (i,val) in self.v_energy.items()])
    self.v_lmap=dict([(imap[i],val) for (i,val) in self.v_lmap.items()])
    self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items()])
    return imap

# ------------------------------------------------------------------------------

  def num_levels(self):
    """ Return number of states. """
    return len(self.v_energy)

# ------------------------------------------------------------------------------

  def nuclear_charge(self):
    """ Return nuclear charge. """
    return self.v_zeff

# ------------------------------------------------------------------------------

  def set_nuclear_charge(self,zeff):
    """ Specify nuclear charge. """
    if zeff is None:
      self.v_zeff=None
    else:
      self.v_zeff=int(zeff)

# ------------------------------------------------------------------------------

  def unset_nuclear_charge(self):
    """ Remove nuclear charge specification. """
    self.v_zeff=None

# ------------------------------------------------------------------------------

  def num_electrons(self):
    """ Return the number of electrons. """
    if self.v_nelec == 0: raise ErrorRAP('no levels present')
    return self.v_nelec

# ------------------------------------------------------------------------------

  def core_energy(self):
    """ Return total energy of all electrons relative to bare nucleus. """
    return self.v_core_energy

# ------------------------------------------------------------------------------

  def set_core_energy(self,core):
    """ Set the core energy in Ry. """
    if core is None:
      self.v_core_energy=None
    else:
      self.v_core_energy=float(core)

# ------------------------------------------------------------------------------

  def unset_core_energy(self):
    """ Remove core energy specification. """
    self.v_core_energy=None

# ------------------------------------------------------------------------------

  def ionization_threshold(self,idx=None):
    """ Return ionization potential of given level (default: ground). """
    if self.v_ion_thresh is None or idx is None: return self.v_ion_thresh
    if idx not in self:
      raise ErrorRAP('level index out of range')
    return self.v_ion_thresh-self.energy(idx)

# ------------------------------------------------------------------------------

  def ionization_potential(self,idx=None):
    """ Return ionization potential of given level (default: ground). """
    return self.ionization_threshold(idx)

# ------------------------------------------------------------------------------

  def set_ionization_threshold(self,ionpot):
    """ Set ionization potential of ground state. """
    if ionpot is None:
      self.v_ion_thresh=None
    else:
      self.v_ion_thresh=float(ionpot)

# ------------------------------------------------------------------------------

  def set_ionization_potential(self,ionpot):
    """ Set ionization potential of ground state. """
    self.set_ionization_threshold(ionpot)

# ------------------------------------------------------------------------------

  def unset_ionization_threshold(self):
    """ Remove ionization potential specification. """
    self.v_ion_thresh=None

# ------------------------------------------------------------------------------

  def parity(self,idx,imap=False):
    """ Return parity of level with given index. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_parity[tidx]

# ------------------------------------------------------------------------------

  def energy_range(self,levlist=None,imap=False):
    """ Return tuple containing the minimum and maximum energies. """
    if imap is True: levlist=map(self._mapme,levlist)
    ens=[self.v_energy[i] for i in levlist if i is not None]
    return (min(ens),max(ens))

# ------------------------------------------------------------------------------

  def energy(self,idx,imap=False):
    """ Return energy of level with given index. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_energy[tidx]

# ------------------------------------------------------------------------------

  def set_energy(self,idx,energy,imap=False):
    """ Set energy of given level to new value. """
    tidx=self._mapme(idx,imap)
    self.v_energy[tidx]=energy

# ------------------------------------------------------------------------------

  def level(self,idx,imap=False):
    """ Return all information for single level (par,E) """
    tidx=self._mapme(idx,imap)
    if tidx is None: return (None,None)
    out=(self.v_parity[tidx], self.v_energy[tidx])
    return out

# ------------------------------------------------------------------------------

  def label(self,idx,imap=False,spacer=None):
    """ Return level label; optional spacer can be used instead of space. """
    raise 'label() is undefined'

# ------------------------------------------------------------------------------

  def weight(self,idx,imap=False):
    """ Return weight of level. """
    raise 'weight() is undefined'

# ------------------------------------------------------------------------------

  def wavelength(self,idx,jdx,imap=False):
    """ Return wavelength in Angstroms calculated from energies. """
    tidx=self._mapme(idx,imap)
    tjdx=self._mapme(jdx,imap)
    if tidx is None or tjdx is None: return None
    ei=self.energy(tidx)
    ej=self.energy(tjdx)
    return HC_RYDANG/abs(ei-ej)

# ------------------------------------------------------------------------------

  def shift_energy(self,en):
    """ Shift all energies by specified amount. """
    ten=float(en)
    for i in self: self.v_energy[i]+=ten

# ------------------------------------------------------------------------------

  def pfilter(self,par,subset=None,imap=False):
    """ Return a list of levels with parity par """
    if par not in [0,1]:
      raise ErrorRAP('invalid parity')
    out=[]
    indices=self.indices()
    if subset is not None: indices=subset
    for i in indices:
      if self.parity(i) == par: out.append(i)
    if imap: out=self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def efilter(self,emin,emax,subset=None,imap=False):
    """ Return a list of levels between emin and emax """
    if emax < emin: emin,emax=emax,emin
    out=[]
    indices=self.indices()
    if subset is not None: indices=subset
    for i in indices:
      e=self.energy(i)
      if e >= emin and e <= emax: out.append(i)
    if imap: out=self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def set_lmap(self,other):
    raise 'set_lmap() undefined'

# ------------------------------------------------------------------------------

  def lmap(self,idx=None):
    """ 
    Return mapped indices.

    If no argument is passed, then return the entire mapping dictionary.  
    Otherwise, the argument can be an integer, list of integers, 2-element
    tuple, or list of 2-element tuples.  The last two items are for mapping
    transitions.
    """
    if idx is None: return self.v_lmap
    if len(self.v_lmap) == 0: return idx
    ttype=type(idx)
    chk=idx
    if ttype == type(0) or ttype == type((1,)): chk=[chk]
    out=[]
    if len(chk) > 0 and isinstance(chk[0],tuple):
      for (i,j) in chk:
        if i in self.v_lmap and j in self.v_lmap:
          out.append((self.v_lmap[i],self.v_lmap[j]))
        else:
          out.append(None)
#      out=[(self.v_lmap[i],self.v_lmap[j]) if i in self.v_lmap and \
#           j in self.v_lmap else None for (i,j) in chk]
    else:
      for i in chk:
        if i in self.v_lmap:
          out.append(self.v_lmap[i])
        else:
          out.append(None)
#      out=[self.v_lmap[i] if i in self.v_lmap else None for i in chk]
    if ttype == type(0) or ttype == type((1,)): out=out[0]
    return out

# ------------------------------------------------------------------------------

  def lmap_reverse(self,idx=None):
    """ 
    Return original indices given mapped values.

    If no argument is passed, then return the entire mapping dictionary.  
    Otherwise, the argument can be an integer, list of integers, 2-element
    tuple, or list of 2-element tuples.  The last two items are for mapping
    transitions.
    """
    if idx is None: return self.v_lmap_r
    if len(self.v_lmap_r) == 0: return idx
    ttype=type(idx)
    chk=idx
    if ttype == type(0) or ttype == type((1,)): chk=[chk]
    out=[]
    if len(chk) > 0 and isinstance(chk[0],tuple):
      for (i,j) in chk:
        if i in self.v_lmap_r and j in self.v_lmap_r:
          out.append((self.v_lmap_r[i],self.v_lmap_r[j]))
        else:
          out.append(None)
#      out=[(self.v_lmap_r[i],self.v_lmap_r[j]) if i in self.v_lmap_r and \
#           j in self.v_lmap_r else None for (i,j) in chk]
    else:
      for i in chk:
        if i in self.v_lmap_r:
          out.append(self.v_lmap_r[i])
        else:
          out.append(None)
#      out=[self.v_lmap_r[i] if i in self.v_lmap_r else None for i in chk]
    if ttype == type(0) or ttype == type((1,)): out=out[0]
    return out

# ------------------------------------------------------------------------------

  def add_level(self):
    raise 'add_level() undefined'

# ------------------------------------------------------------------------------

  def add_levels(self):
    raise 'add_level() undefined'

################################################################################
class rap_trans:
  """ 
  Class for storing transition data; also serves as a superclass for
  specific transition data, e.g. A-values.
  """
  def __init__(self,levs=None):
    self.clear()
    if levs is not None: self.set_levs(levs)

# ------------------------------------------------------------------------------

  def __len__(self):
    return self.num_transitions()

# ------------------------------------------------------------------------------

  def __getitem__(self,key):
    return self.v_data.__getitem__(key)

# ------------------------------------------------------------------------------

  def __contains__(self,key):
    return self.v_data.__contains__(key)

# ------------------------------------------------------------------------------

  def __iter__(self):
    return self.v_data.__iter__()

# ------------------------------------------------------------------------------

  def _relimit(self):
    """ Remove transitions not in given initial/final state limits """
    if self.v_iminlev is not None:
      for (i,j) in self.v_data.keys():
        if i < self.v_iminlev: del self.v_data[(i,j)]
    if self.v_imaxlev is not None:
      for (i,j) in self.v_data.keys():
        if i > self.v_imaxlev: del self.v_data[(i,j)]
    if self.v_fminlev is not None:
      for (i,j) in self.v_data.keys():
        if j < self.v_fminlev: del self.v_data[(i,j)]
    if self.v_fmaxlev is not None:
      for (i,j) in self.v_data.keys():
        if j > self.v_fmaxlev: del self.v_data[(i,j)]

# ------------------------------------------------------------------------------

  def clear(self,keeplevs=False):
    """ Reset data. """
    if not keeplevs: self.v_levs=None
    self.v_iminlev=None     # lower limit for initial state index
    self.v_imaxlev=None     # upper limit for initial state index
    self.v_fminlev=None     # lower limit for final state index
    self.v_fmaxlev=None     # upper limit for final state index
    self.v_data={}

# ------------------------------------------------------------------------------

  def reindex(self,remap):
    """ 
    Reassign level indices.  Indices are remapped according to dictionary,
    remap, old (key) -> new (value).  
    """
    indices=[]
    for (i,j) in self:
      if i not in indices: indices.append(i)
      if j not in indices: indices.append(j)
    err=False
    for i in indices:
      if i not in remap:
        err=True
        break
    if err:
      raise ErrorRAP("all indices must be represented using reindex")
    tdat={}
    for (i,j),val in self.v_data.iteritems(): tdat[(remap[i],remap[j])]=val
    self.v_data=tdat
    if self.v_iminlev is not None: self.v_iminlev=remap[self.v_iminlev]
    if self.v_imaxlev is not None: self.v_imaxlev=remap[self.v_imaxlev]
    if self.v_fminlev is not None: self.v_fminlev=remap[self.v_fminlev]
    if self.v_fmaxlev is not None: self.v_fmaxlev=remap[self.v_fmaxlev]
    return remap

# ------------------------------------------------------------------------------

  def set_levs(self,levs):
    """ Set rap_states object. """
    if not isinstance(levs,rap_states):
      msg='levs must be a rap_states instance (or child instance)'
      raise ErrorRAP(msg)
    self.v_levs=levs

# ------------------------------------------------------------------------------

  def levs(self):
    """ Return rap_states object. """
    return self.v_levs

# ------------------------------------------------------------------------------

  def num_levels(self):
    """ Return number of levels. """
    return len(self.v_levs)

# ------------------------------------------------------------------------------

  def num_transitions(self):
    """ Return number of transitions. """
    return len(self.v_data)

# ------------------------------------------------------------------------------

  def initial_limit(self,imap=False):
    """ Return allowed range for inital state indices; None for no limit. """
    imin=self.v_iminlev
    imax=self.v_imaxlev
    if imap is False: return (imin,imax)
    if imin is None and imax is None: return (imin,imax)
    tmin=1
    if imin is not None: tmin=imin
    tmax=self.num_levels()
    if tmax is not None: tmax=imax
    levlist=[self.v_levs.lmap_reverse(i) for i in range(tmin,tmax+1)]
    tmin=min(levlist)
    tmax=max(levlist)
    if imin is None: tmin=None
    if imax is None: tmax=None
    return (tmin,tmax)

# ------------------------------------------------------------------------------

  def final_limit(self,imap=False):
    """ Return allowed range for final state indices; None for no limit. """
    fmin=self.v_fminlev
    fmax=self.v_fmaxlev
    if imap is False: return (fmin,fmax)
    if fmin is None and fmax is None: return (fmin,fmax)
    tmin=1
    if fmin is not None: tmin=fmin
    tmax=self.num_levels()
    if tmax is not None: tmax=fmax
    levlist=[self.v_levs.lmap_reverse(i) for i in range(tmin,tmax+1)]
    tmin=min(levlist)
    tmax=max(levlist)
    if fmin is None: tmin=None
    if fmax is None: tmax=None
    return (tmin,tmax)

# ------------------------------------------------------------------------------

  def set_initial_limit(self,minlev=None,maxlev=None,imap=False):
    """ Set allowed range for inital state indices; None for no limit. """
    tmin=None
    if minlev is not None: tmin=self.v_levs._mapme(minlev,imap)
    tmax=None
    if maxlev is not None: tmax=self.v_levs._mapme(maxlev,imap)
    if tmin is not None and tmax is not None:
      if tmin > tmax: tmin,tmax=tmax,tmin
    indices=self.v_levs.indices()
    if tmin is not None and tmin < min(indices):
      msg='minlev out of range'
      raise ErrorRAP(msg)
    if tmax is not None and tmax > max(indices):
      msg='maxlev out of range'
      raise ErrorRAP(msg)
    self.v_iminlev=tmin
    self.v_imaxlev=tmax
    self._relimit()

# ------------------------------------------------------------------------------

  def set_final_limit(self,minlev=None,maxlev=None,imap=False):
    """ Set allowed range for final state indices; None for no limit. """
    tmin=None
    if minlev is not None: tmin=self.v_levs._mapme(minlev,imap)
    tmax=None
    if maxlev is not None: tmax=self.v_levs._mapme(maxlev,imap)
    if tmin is not None and tmax is not None:
      if tmin > tmax: tmin,tmax=tmax,tmin
    indices=self.v_levs.indices()
    if tmin is not None and tmin < min(indices):
      msg='minlev out of range'
      raise ErrorRAP(msg)
    if tmax is not None and tmax > max(indices):
      msg='maxlev out of range'
      raise ErrorRAP(msg)
    self.v_fminlev=tmin
    self.v_fmaxlev=tmax
    self._relimit()

# ------------------------------------------------------------------------------

  def initial_indices(self,imap=False):
    """ Return list of active indices. """
    out=list(set([sublist[0] for sublist in self.v_data.keys()]))
    (imin,imax)=self.initial_limit()
    if imin is not None: out=filter(lambda x: x >= imin,out)
    if imax is not None: out=filter(lambda x: x <= imax,out)
    if imap: out=self.v_levs.lmap_reverse(out)
    return sorted(out)

# ------------------------------------------------------------------------------

  def final_indices(self,imap=False):
    """ Return list of active indices. """
    out=list(set([sublist[1] for sublist in self.v_data.keys()]))
    (imin,imax)=self.initial_limit()
    if imin is not None: out=filter(lambda x: x >= imin,out)
    if imax is not None: out=filter(lambda x: x <= imax,out)
    if imap: out=self.v_levs.lmap_reverse(out)
    return sorted(out)

# ------------------------------------------------------------------------------

  def indices(self,imap=False):
    """ Return list of active indices. """
    out=self.initial_indices()
    out.extend(self.final_indices())
    out=list(set(out))
    if imap: out=self.v_levs.lmap_reverse(out)
    return sorted(out)

# ------------------------------------------------------------------------------

  def initial_filter(self,istates,imap=False):
    """ Return list of active indices. """
    tmp=istates
    if not isinstance(tmp,list): tmp=[tmp]
    out=[key for key in self.v_data.keys() if key[0] in tmp]
    if imap: out=self.v_levs.lmap_reverse(out)
    return sorted(out)

# ------------------------------------------------------------------------------

  def final_filter(self,istates,imap=False):
    """ Return list of active indices. """
    tmp=istates
    if not isinstance(tmp,list): tmp=[tmp]
    out=[key for key in self.v_data.keys() if key[1] in tmp]
    if imap: out=self.v_levs.lmap_reverse(out)
    return sorted(out)

# ------------------------------------------------------------------------------

  def transitions(self,imap=False):
    """ Return list of transition tuples: (lower,upper). """
    out=self.v_data.keys()
    if imap: out=self.v_levs.lmap_reverse(out)
    return sorted(out)

# ------------------------------------------------------------------------------

  def data(self,i,j=None,imap=False):
    """ 
    Return the data value for the i-j transition.
    
    Note: it is assumed that state i has a lower energy than state j.
    """
    ii,jj=i,j
    if jj is None and isinstance(ii,tuple): (ii,jj)=ii
    trans=(ii,jj)
    if imap:
      trans=self.v_levs.lmap(trans)
      if trans is None: raise ErrorRAP('invalid mapping')
    if trans not in self.v_data: 
      (ii,jj)=trans
      raise ErrorRAP('Data value for %i-%i transition not found.'%(ii,jj))
    return self.v_data[trans]

# ------------------------------------------------------------------------------

  def set_data(self,i,j=None,value=None,imap=False):
    """ 
    Set data value for the i-j transition.
    
    Note: it is assumed that state i has a lower energy than state j.
    """
    if isinstance(i,tuple):
      (ii,jj)=i
      tval=j
    else:
      ii,jj=i,j
      tval=value
    if tval is None:
      raise ErrorRAP("need to give a new value")
    if ii not in self.levs():
      raise ErrorRAP("Invalid state index given: %s"%ii)
    if jj not in self.levs():
      raise ErrorRAP("Invalid state index given: %s"%jj)
    trans=(ii,jj)
    if imap:
      trans=self.v_levs.lmap(trans)
      if trans is None: raise ErrorRAP('invalid mapping')
    self.v_data[trans]=tval

# ------------------------------------------------------------------------------

  def wavelength(self,i,j,imap=False):
    """ Return wavelength calculated from level structure in Angstroms. """
    return self.levs().wavelength(i,j,imap=imap)

################################################################################

################################################################################
class rap_states_ls(rap_states):
  """ Store set of LS states. """
  def __init__(self,filename=None,**args):
    self.clear()
    if filename: self.loadfile(filename,**args)

# ------------------------------------------------------------------------------

  def _reprocess(self):
    """ Redefine senority-like tags. """
    if len(self) == 0: return
    indices=zip(*sorted(zip(self.v_energy.values(),self.v_energy.keys())))[1]
    allowedterms={}

    have={}
    for idx in indices:
      (conf,term,sen,par,en)=self.level(idx)

      # check for legal term
      if conf not in allowedterms: allowedterms[conf]=conf.allowed_terms()
      if term not in allowedterms[conf]:
        raise ErrorRAP('invalid term for given configuration')

      # get and store seniority
      key=(conf,term)
      if key not in have: have[key]=0
      have[key]+=1
      self.v_seniority[idx]=have[key]
      if have[key] > allowedterms[conf][term]:
        raise ErrorRAP('too many terms in configuration')

# ------------------------------------------------------------------------------

  def clear(self):
    """ Erase all data. """
    rap_states.clear(self)
    self.v_cup='LS'
    self.v_conf={}
    self.v_term={}
    self.v_seniority={}

# ------------------------------------------------------------------------------

  def clear_levels(self):
    """ Just erase level data. """
    rap_states.clear_levels(self)
    self.v_cup='LS'
    self.v_conf={}
    self.v_term={}
    self.v_seniority={}

# ------------------------------------------------------------------------------

  def loadfile(self,filename,**args):
    """ Load in level information from level data file. """
    ft=rlf_filetype(filename)
    if ft == 'adf04':
      self.load_adf04(filename,**args)
    elif ft == 'olg':
      self.load_olg(filename,**args)
    elif ft == 'x006':
      self.load_xstar(filename,**args)
    elif ft == 'tiptop':
      self.load_tiptop(filename,**args)
    else:
      raise ErrorRAP('do not know how to read file type')

# ------------------------------------------------------------------------------

  def reindex(self,remap=None):
    """ 
    Reassign level indices.  Indices are remapped according to dictionary,
    remap, old (key) -> new (value).  If remap is None, then reindex according
    to energy order.
    """
    imap=rap_states.reindex(self,remap)
    self.v_conf=dict([(imap[i],val) for (i,val) in self.v_conf.items()])
    self.v_term=dict([(imap[i],val) for (i,val) in self.v_term.items()])
    self.v_seniority=dict([(imap[i],val) for (i,val) in \
                          self.v_seniority.items()])

# ------------------------------------------------------------------------------

  def configurations(self,levlist=None,imap=False):
    """ Return list of all configurations. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self._mapme,levlist)
    out=[]
    for i in levlist:
      if self.v_conf[i] in out: continue
      out.append(self.v_conf[i])
    return out

# ------------------------------------------------------------------------------

  def configuration(self,idx,imap=False):
    """ Return configuration of level with given index. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_conf[tidx]

# ------------------------------------------------------------------------------

  def terms(self,levlist=None,imap=False):
    """ Return list of all terms. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self.mapme,levlist)
    out=[]
    for i in levlist:
      if self.v_term[i] in out: continue
      out.append(self.v_term[i])
    return out

# ------------------------------------------------------------------------------

  def term(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_term[tidx]

# ------------------------------------------------------------------------------

  def seniority(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_seniority[tidx]

# ------------------------------------------------------------------------------

  def level(self,idx,imap=False):
    """ Return all information for single level (conf,term,alpha,par,E) """
    tidx=self._mapme(idx,imap)
    if tidx is None: return (None,None,None,None,None)
    out=(self.v_conf[tidx],self.v_term[tidx],self.v_seniority[tidx],\
         self.v_parity[tidx], self.v_energy[tidx])
    return out

# ------------------------------------------------------------------------------

  def label(self,idx,imap=False,spacer=None):
    """ Return level label; optional spacer can be used instead of space. """
    tidx=self._mapme(idx,imap)
    out=str(self.v_conf[tidx])+' '+self.v_term[tidx]
    if self.v_seniority[tidx] > 1: out+='#'+str(self.v_seniority[tidx])
    if spacer is not None: out=out.replace(' ',spacer)
    return out

# ------------------------------------------------------------------------------

  def weight(self,idx,imap=False):
    """ Return weight of level. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return rlf_termweight(self.v_term[tidx])

# ------------------------------------------------------------------------------

  def add_level(self,conf,term,energy,idx=False):
    """ Add new level to list """
    if idx is not False and idx in self:
      raise ErrorRAP('index already exists')
    tidx=idx
    if tidx is False: 
      tidx=1
      if len(self) > 0: tidx=max(self.indices())+1
    tconf=conf
    if not isinstance(tconf,rap_config): tconf=rap_config(tconf)
    if self.num_electrons() is None:
      self.v_nelec=tconf.num_electrons()
    else:
      if tconf.num_electrons() != self.num_electrons():
        raise ErrorRAP('inconsistent number of electrons')
    (g,l)=rlf_termdecode(term)     # will raise error if bad term string
    self.v_conf[tidx]=tconf
    self.v_term[tidx]=term
    self.v_seniority[tidx]=-1
    self.v_parity[tidx]=tconf.parity()
    self.v_energy[tidx]=energy
    self._reprocess()

# ------------------------------------------------------------------------------

  def add_levels(self,toadd,idxstart=None):
    """ 
    Add set of levels to structure. 

    The input list, toadd, contains tuples with: energy (Ry), index,
    configuration, and term.  If the indexing is to be automated, then set 
    index in the tuples to -1.
    """
    if len(toadd) == 0: return
    toadd.sort()
    idxcur=1
    if len(self) > 0: idxcur=max(self.indices())+1
    if idxstart is not None:
      if idxstart < idxcur: raise ErrorRAP('invalid starting index')
      idxcur=idxstart
    tmp={}
    nechk=None
    for (en,idx,conf,term) in toadd:
      tidx=idx
      if tidx in self: raise ErrorRAP('duplicate index')
      if tidx in tmp: raise ErrorRAP('duplicate index')
      if tidx < 0:
        tidx=idxcur
        idxcur+=1
      tconf=conf
      if not isinstance(conf,rap_config): tconf=rap_config(conf)
      if nechk is None: nechk=tconf.num_electrons()
      if tconf.num_electrons() != nechk:
        raise ErrorRAP('inconsistent number of electrons')
      (g,l)=rlf_termdecode(term)     # will raise error if bad term string
      tmp[tidx]=(en,tconf,term)
    self.v_nelec=nechk
    for idx,(en,conf,term) in tmp.iteritems():
      self.v_conf[idx]=conf
      self.v_term[idx]=term
      self.v_seniority[idx]=-1
      self.v_parity[idx]=conf.parity()
      self.v_energy[idx]=en
    self._reprocess()

# ------------------------------------------------------------------------------

  def set_lmap(self,other,simple=False):
    """ 
    Map to other rap_states_ls instance. 
    
    By default, energies are used when differing numbers of duplicate terms
    for a given configuration are present (one set of energies being 
    incomplete).  The keyword, simple, can be set to True, however, to use
    a simple mapping without using energies where terms are associated by
    order of appearance.
    """
    if not isinstance(other,rap_states_ls):
      raise ErrorRAP('argument must be of type rap_states_ls')
    if self.num_electrons() != other.num_electrons():
      raise ErrorRAP('inconsistent number of electrons')

    if simple:
      tmap={}
      for i in other:
        k=None
        for j in self:
          if self[j] == other[i]:
            k=j
            break
        tmap[i]=k
      self.v_lmap=tmap
      self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items()])
      return

    # group level indices by (conf,term); seniority numbers not trusted
    fromlist={}
    for i in self:
      key=(self.configuration(i),self.term(i))
      if key not in fromlist: fromlist[key]=[]
      fromlist[key].append(i)
    tolist={}
    for j in other:
      key=(other.configuration(j),other.term(j))
      if key not in tolist: tolist[key]=[]
      tolist[key].append(j)

    self.v_lmap=self._map_by_energy(other,fromlist,tolist)
    for j in other:
      if j not in self.v_lmap: self.v_lmap[j]=None
    self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items() if v is not None])
    for i in self:
      if i not in self.v_lmap_r: self.v_lmap_r[i]=None

# ------------------------------------------------------------------------------

  def cfilter(self,clist,subset=None,imap=False):
    """ Filter list of states based on configurations in clist. """
    if not isinstance(clist,list) and not isinstance(clist,tuple): clist=[clist]
    tclist=[]
    for c in clist:
      tc=c
      if not isinstance(tc,rap_config): tc=rap_config(tc)
      tclist.append(tc)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.configuration(i) in tclist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def tfilter(self,tlist,alpha=None,subset=None,imap=False):
    """ 
    Filter list of states based on terms in tlist; allow for optional input
    for the seniority (alpha).
    """
    if not isinstance(tlist,list) and not isinstance(tlist,tuple): tlist=[tlist]
    if alpha is not None:
      if not isinstance(alpha,list) and not isinstance(alpha,tuple): 
        alpha=[alpha]

    ttlist=[]
    for t in tlist:
      tt=t
      if isinstance(tt,tuple): tt=rlf_termencode(tt[0],tt[1])
      ttlist.append(tt)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if alpha is not None and self.seniority(i) not in alpha: continue
      if self.term(i) in ttlist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def symfilter(self,g,l=None,p=None,subset=None,imap=False):
    """ Filter list of states based on symmetry. """
    if l is None and p is None:
      try:
        (gg,ll,pp)=g
      except:
        msg='do not understand symmetry'
        raise ErrorRAP(msg)
    else:
      gg=g
      ll=l
      pp=p

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      (tg,tl)=rlf_termdecode(self.term(i))
      if tg != gg or tl != ll: continue
      tp=self.parity(i)
      if tp != pp: continue
      out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def LSJ_split(self):
    """ 
    Split LS energy level data in LSJ states.  Energies for each level are the
    same as the original term, but the order is determined by Hund's Rule #3.
    """
    toadd=[]
    idx=0
    for i in self:
      (conf,term,alfa,p,en)=self.level(i)

      # determine order of levels
      ascend=False
      for ss in conf:
        if conf[ss] < rlf_ssoccupancy(ss)/2: 
          ascend=True

      j2vals=rlf_term_2jvals(term)
      if not ascend: j2vals.reverse()

      for j2 in j2vals: 
        idx+=1
        toadd.append((en,idx,conf,term,j2))

    # assemble LSJ object
    out=rap_states_lsj()
    out.set_nuclear_charge(self.nuclear_charge())
    out.set_core_energy(self.core_energy())
    out.set_ionization_threshold(self.ionization_threshold())
    out.add_levels(toadd)
    return out

# ------------------------------------------------------------------------------

  def load_olg(self,filename):
    """ load information from olg file """
    fin=open(filename,'r')
    self.v_levels=[]

    # first get closed subshells and number of configurations
    closed=[]
    while 1:
      line=fin.readline()
      tmp=re.findall(r'IN(.+)CONFIGURATIONS',line)
      if line == '': RuntimeError('configuration list not found')
      if tmp: break
    nconf=int(tmp[0])
    while 1:
      tmp=re.findall(r'(\d+)\s+(\d+)',line)
      if len(tmp) == 0: break
      closed.extend(tmp)
      line=fin.readline()

    # get list of all open subshells
    while 1:
      line=fin.readline()
      if line == '': raise ErrorRAP('could not find sub-shell list')
      if 'INTERNAL ORBITALS' in line: break
    torbs=[]
    while 1:
      line=fin.readline()
      if 'CONFIGURATION' in line: break
      if 'N L' in line: 
        tmp=re.findall(r'(\d+)\s+(\d+)',line)
        if tmp: torbs.extend(tmp)
    orbs=[]
    for o in torbs:
      if o not in closed: orbs.append(o)
    base={}
    for o in orbs: base[o]=0

    ne=None
    conflist=[]
    for iconf in range(nconf):
      while 1:
        if line == '': return False
        tmp=re.findall(r'COMBINATIONS(.+)',line)
        if tmp: break
        line=fin.readline()
      tmp=tmp[0]
      while 1:
        line=fin.readline()
        if 'SLATER' in line: break
        tmp+=line[:-1]
      tconf=base.copy()
      nl=re.findall(r'(\d+)\s+(\d+)',tmp)
      if not nl: raise ErrorRAP('bad sub-shell')
      for ss in nl: tconf[ss]+=1
      aconf=''
      for o in orbs:
        ss=str(o[0])+SPECLETTERS[int(o[1])].lower()
        n=tconf[o]
        aconf+=ss+str(n)+' '
      conflist.append(rap_config(aconf,ne=ne))
      if ne is None: ne=conflist[0].num_electrons()

    # get nuclear charge
    while 1:
      if line == '':
        raise ErrorRAP('premature end of file')
      line=fin.readline()
      if 'ATOMIC NUMBER' in line: break
    tmp=re.findall(r'ATOMIC NUMBER\s*(\d+)',line)
    zeff=int(tmp[0])

    # get level information
    while 1:
      line=fin.readline()
      if line == '': raise ErrorRAP('no level data found')
      if '(EI-E1)/RY' in line: break
    tmp=re.findall(r'(-?\d+\.\d+)\s*$',line)
    core_energy=None
    if tmp: core_energy=float(tmp[0])
    toadd=[]
    while 1:
      line=fin.readline()
      tmp=line.split()
      if tmp == []: break
      try:
        idx=int(tmp[0])
      except:
        break
      term=rlf_termencode(abs(int(tmp[3])),int(tmp[4]))
      conf=conflist[int(tmp[5])-1]
      energy=float(tmp[7])
      toadd.append((energy,idx,conf,term))
    fin.close()
    self.set_nuclear_charge(zeff)
    self.set_core_energy(core_energy)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_adf04(self,filename):
    """ Load in level information from an ADAS adf04 file """
    fin=open(filename,'r')

    # basic info
    line=fin.readline()
    tmp=re.findall(r'([\d.]+)',line)
    z=int(tmp[1])
    q=int(tmp[2])
    ip=float(tmp[3])
    ne=z-q+1

    # read in level data
    isic=True
    isls=True
    isca=True
    toadd=[]
    while 1:
      line=fin.readline()
      if '-1' in line: break
      tmp=re.findall(r'(\d+)(.+)\((.+)\)(.+)\((.+)\)\s+(\d+[\.\d]*)',line)[0]
      idx=int(tmp[0])
      conf=rap_config(tmp[1],ne=ne)
      g=int(tmp[2])
      l=int(tmp[3])
      j2=int(2.*float(tmp[4]))
      e=float(tmp[5])/INVERSECM
      toadd.append((e,idx,conf,rlf_termencode(g,l)))

      cwgt=conf.weight()
      lwgt=g*(2*l+1)
      wgt=j2+1
      if wgt != cwgt: isca=False
      if wgt != lwgt: isls=False
      if wgt not in  range(abs(2*l-g+1)+1,2*l+g+1,2): isic=False
    fin.close()

    # check for correct coupling scheme
    if not isls:
      if isca: 
        raise ErrorRAP('adf04 file seems to be configuration averaged')
      elif isic:
        raise ErrorRAP('adf04 file seems to be J-resolved')
      else:
        raise ErrorRAP('adf04 file not LS or a mistake in the file')

    self.set_nuclear_charge(z)
    self.set_ionization_threshold(ip)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_tiptop(self,filename,reorder=False):
    """ 
    Read energy levels from TIPTOP energy file.

    reorder - If True, reindex indices according to energy order.  
    """
    fin=open(filename,'r')

    # read header and get column meanings
    line=fin.readline()
    hline=fin.readline()
    line=fin.readline()

    tmp=hline.split()
    ncol=len(tmp)
    id_i=tmp.index('i')
    id_nz=tmp.index('NZ')
    id_ne=tmp.index('NE')
    id_islp=tmp.index('iSLP')
    id_eryd=tmp.index('E(RYD)')
    id_iconf=tmp.index('iCONF')

    # There can be spaces in the configuration strings, which can
    # mess up indices for columns after the configuration and make
    # getting the configuration string complicated; thus the following
    id_iconf_end=id_iconf-ncol
    if id_i > id_iconf: id_i=id_i-ncol
    if id_nz > id_iconf: id_nz=id_nz-ncol
    if id_ne > id_iconf: id_ne=id_ne-ncol
    if id_islp > id_iconf: id_islp=id_islp-ncol
    if id_eryd > id_iconf: id_eryd=id_eryd-ncol
 
    # read in data
    nz=None
    ip=None
    toadd=[]
    while 1:
      line=fin.readline()
      if '=====' in line: break

      tmp=line.split()
      idx=int(tmp[id_i])
      tnz=int(tmp[id_nz])
      islp=tmp[id_islp]
      eryd=float(tmp[id_eryd])
      confstr=' '.join(tmp[id_iconf:id_iconf_end+1])
      conf=rap_config(confstr)
      if idx == 1: ip=abs(eryd)
      if reorder: idx=None

      g=islp[0]
      l=SPECLETTERS[int(islp[1])]
      term=str(g)+l

      if nz is None: nz=tnz
      toadd.append((eryd,idx,conf,term))
    fin.close()

    # ground state energy is given wrt ionization threshold; need to shift
    newadd=[]
    for (eryd,idx,conf,term) in toadd:
      enew=eryd+ip
      if idx == 1: enew=0.0       # be sure to avoid precision errors
      newadd.append((enew,idx,conf,term))

    self.set_ionization_threshold(ip)
    self.set_nuclear_charge(nz)
    self.add_levels(newadd)

# ------------------------------------------------------------------------------

  def load_xstar(self,filename):
    """ 
    Load level information from xstar type 006 file.
    """
    # column id's
    id_energy=7
    id_j2=8
    id_ip=10
    id_g=12
    id_l=13
    id_z=14
    id_ion=16
    id_id=15
    id_conf=17

    # read level list from file
    ion=None
    ne=None
    ldata=[]
    have=[]
    xid0=None
    fin=open(filename,'r')
    for line in fin:
      if line.strip() == '': continue
      data=line.split()

      confstr=data[id_conf-1]
      xid=int(data[id_id-1])
      ion=int(data[id_ion-1])
      g=abs(int(data[id_g-1]))
      l=int(data[id_l-1])
      j2=int(float(data[id_j2-1]))-1
      w=g*(2*l+1)
      if g == 1 or l == 0:
        cup='XX'
      elif w < j2+1:
        cup='CA'
        if 'CA' not in have: have.append('CA')
      elif w == j2+1:
        cup='LS'
        if 'LS' not in have: have.append('LS')
      else:
        j2min=abs(g-1-2*l)
        j2max=g-1+2*l
        okay=False
        for j2t in range(j2min,j2max+1,2):
          if j2t == j2: okay=True
        if okay:
          cup='IC'
          if 'IC' not in have: have.append('IC')
        else:
          cup='??'
      e=float(data[id_energy-1])
      if e == 0.0:
        ne=rap_config(confstr).num_electrons()
        zeff=int(data[id_z-1])
        ion_thresh=float(data[id_ip-1])
      chk=confstr.lower()
      if 'su' in chk or 'sp' in chk: continue 
      if 'co' in chk or 'cn' in chk: continue
      if 'fake' in chk: continue
      ldata.append([e/RYDBERG,cup,confstr,g,l,j2,xid])
    fin.close()
    if ne is None: raise RuntimeError('cannot find ground state')

    # id's for ldata
    id_e=0
    id_cup=1
    id_conf=2
    id_g=3
    id_l=4
    id_j2=5
    id_xid=6

    # get number of electrons from ground state and standardize configurations
    for i in range(len(ldata)):
      ldata[i][id_conf]=rap_config(ldata[i][id_conf],ne=ne)

    # get total weight of each configuration
    weight={}
    wtot={}
    cups={}
    for i in range(len(ldata)):
      conf=str(ldata[i][id_conf])
      if conf not in wtot: wtot[conf]=ldata[i][id_conf].weight()
      if conf not in cups: cups[conf]=[]
      if ldata[i][id_cup] not in cups[conf]: cups[conf].append(ldata[i][id_cup])
      if ldata[i][id_cup] in ['XX','LS']:
        if conf not in weight:
          weight[conf]=ldata[i][id_g]*(2*ldata[i][id_l]+1)
        else:
          weight[conf]+=ldata[i][id_g]*(2*ldata[i][id_l]+1)
      else:
        if conf not in weight:
          weight[conf]=(ldata[i][id_j2]+1)
        else:
          weight[conf]+=(ldata[i][id_j2]+1)

    # check if configurations have a common coupling scheme
    for c in cups.keys():
      tmp=cups[c][:]
      if 'XX' in tmp: tmp.remove('XX')
      if len(tmp) > 1:
        raise RuntimeError('mixed coupling scheme not supported')
      if tmp == [] and len(have) == 1: tmp=have[:]
      cups[c]=tmp[:]
      if cups[c] == []: cups[c]=['LS']

    # check that 'LS' configurations have all possible terms
    for c in weight.keys():
      if 'LS' not in cups[c]: continue
      elif weight[c] > wtot[c]:
        raise RuntimeError('configuration, '+str(c)+', has too many terms')

    # add levels having LS coupling
    toadd=[]
    for lev in ldata:
      idx=lev[id_xid]
      conf=lev[id_conf]
      sc=str(conf)
      if 'LS' not in cups[sc] and cups[sc] != ['XX']: continue
      term=rlf_termencode(lev[id_g],lev[id_l])
      energy=lev[id_e]
      toadd.append((energy,idx,conf,term))
    self.set_nuclear_charge(zeff)
    self.set_ionization_threshold(ion_thresh/RYDBERG)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def write_xstar(self,filename,subset=None,precision=5):
    """ 
    Write xstar type 6 file. Must have ionization potential and nuclear charge
    specified.
    """
    if self.ionization_threshold() == None:
      raise ErrorRAP('must define ionization threshold')
    if self.nuclear_charge() == None:
      raise ErrorRAP('must define nuclear charge')
    indices=self.indices()
    if subset is not None: indices=subset

    zeff=self.nuclear_charge()
    ne=self.num_electrons()
    ionno=xstar_ionencode(zeff,ne)

    aform='%'+str(8+precision)+'.'+str(precision)+'E'
    iform='%6i'
    prefix=iform%6+iform%13+iform%0+iform%4+iform%6
    fout=open(filename,'w')
    for i in indices:
      conf=self.configuration(i)
      e=self.energy(i)*RYDBERG
      wgt=float(self.weight(i))
      ip=self.ionization_threshold()*RYDBERG
      np=-1.0
      if ip > 0.0: 
        z1=float(zeff-ne+1)
        de=abs(self.ionization_threshold()-self.energy(i))
        np=z1/sqrt(de)
      nval=0
      if conf is not None: nval=rlf_ssdecode(conf.valence_subshell())[0]
      (g,l)=rlf_termdecode(self.term(i))
      if g is None: 
        g,l=0,-1
      else:
        par=conf.parity()
        if par == 1: g=-g
      label=self.label(i,spacer='.')
      llen=len(label)
      line=prefix+iform%llen+' '+aform%e+aform%wgt+aform%np+aform%ip+' '+\
                  iform%nval+iform%g+iform%l+iform%zeff+iform%i+\
                  iform%ionno+' '+label+' %'
      fout.write(line+'\n')
    fout.close()

################################################################################
class rap_rad(rap_trans):
  """ 
  Class for storing Einstein A coefficients.
  """
  def __init__(self,filename=None,levs=None,**args):
    self.clear()
    if levs is not None: self.set_levs(levs)
    if filename is not None: self.loadfile(filename,**args)

# ------------------------------------------------------------------------------

  def clear(self,keeplevs=False):
    """ Erase all data. """
    rap_trans.clear(self,keeplevs=keeplevs)
    self.v_lifetime={}

# ------------------------------------------------------------------------------

  def loadfile(self,filename,cup=None,iminlev=None,imaxlev=None,\
                                      fminlev=None,fmaxlev=None):
    """ Load radiative data from file. """
    ftype=rlf_filetype(filename)
    if ftype == 'olg':
      self.read_olg(filename,cup=cup,iminlev=iminlev,imaxlev=imaxlev,\
                                     fminlev=fminlev,fmaxlev=fmaxlev)
    elif ftype == 'adf04':
      self.read_adf04(filename,cup=cup,iminlev=iminlev,imaxlev=imaxlev,\
                                       fminlev=fminlev,fmaxlev=fmaxlev)
    elif ftype == 'chianti_wgfa':
      self.read_chianti(filename,iminlev=iminlev,imaxlev=imaxlev,\
                                 fminlev=fminlev,fmaxlev=fmaxlev)
    elif ftype == 'x050':
      self.read_xstar(filename,iminlev=iminlev,imaxlev=imaxlev,\
                                fminlev=fminlev,fmaxlev=fmaxlev)
    elif ftype == 'nist_lines':
      self.read_nist(filename,iminlev=iminlev,imaxlev=imaxlev,\
                              fminlev=fminlev,fmaxlev=fmaxlev)
    elif ftype == 'tiptop_lines':
      self.read_tiptop(filename,iminlev=iminlev,imaxlev=imaxlev,\
                                fminlev=fminlev,fmaxlev=fmaxlev)
    else:
      raise ErrorRAP('unsupported file type')

# ------------------------------------------------------------------------------

  def reindex(self,remap):
    """ 
    Reassign level indices.  Indices are remapped according to dictionary,
    remap, old (key) -> new (value).  
    """
    imap=rap_trans.reindex(self,remap)
    self.v_lifetime=dict([(imap[i],val) for (i,val) in self.v_lifetime.items()])
    return imap

# ------------------------------------------------------------------------------

  def avalue(self,i,j=None,imap=False):
    """ 
    Return the A-value for the i-j transition.
    
    Note: it is assumed that state i has a lower energy than state j.
    """
    return self.data(i,j,imap)

# ------------------------------------------------------------------------------

  def set_avalue(self,i,j,avalue,imap=False):
    """ Add or update an A-value (also updates lifetime). """
    ii,jj=i,j
    if imap: (ii,jj)=self.v_levs.lmap((ii,jj))
    if self.v_levs.energy(ii) > self.v_levs.energy(jj): ii,jj=jj,ii
    (iimin,iimax)=self.initial_limit()
    if ii < iimin and ii > iimax:
      msg='initial index outside of user-defined range: %s'%ii
      raise ErrorRAP(msg)
    (jjmin,jjmax)=self.final_limit()
    if jj < jjmin and jj > jjmax:
      msg='final index outside of user-defined range: %s'%jj
      raise ErrorRAP(msg)
    if (ii,jj) in self:
      tau=self.v_lifetime[jj]
      taup=1./(1./tau-self[(ii,jj)]+avalue)
    else:
      taup=1./avalue
    self.v_data[(ii,jj)]=avalue
    self.v_lifetime[jj]=taup

# ------------------------------------------------------------------------------

  def set_avalues(self,avalues,imap=False):
    """ 
    Efficiently add many new A-values.  Input is a dictionary with keys being
    the transition tuple by index (lower,upper) and the values are the A-values.
    """
    newsum={}
    oldsum={}
    for (i,j),av in avalues.iteritems():
      ii,jj=i,j
      if imap: (ii,jj)=self.v_levs.lmap((ii,jj))
      if self.v_levs.energy(ii) > self.v_levs.energy(jj): ii,jj=jj,ii
      (iimin,iimax)=self.initial_limit()
      if ii < iimin and ii > iimax:
        msg='initial index outside of user-defined range: %s'%ii
        raise ErrorRAP(msg)
      (jjmin,jjmax)=self.final_limit()
      if jj < jjmin and jj > jjmax:
        msg='final index outside of user-defined range: %s'%jj
        raise ErrorRAP(msg)
      if jj not in newsum: 
        newsum[jj]=0.0
        oldsum[jj]=0.0
      newsum[jj]+=av
      if (ii,jj) in self: oldsum[jj]+=self[(ii,jj)]
      self.v_data[(ii,jj)]=av
    for jj in newsum:
      if jj not in self.v_lifetime:
        self.v_lifetime[jj]=1./newsum[jj]
      else:
        self.v_lifetime[jj]=1./(1./self.v_lifetime[jj]-oldsum[jj]+newsum[jj])


# ------------------------------------------------------------------------------

  def gavalue(self,i,j=None,imap=False):
    """ 
    Return the gA-value for the i-j transition.
    
    Note: it is assumed that state i has a lower energy than state j.
    """
    ii,jj=i,j
    if jj is None and isinstance(ii,tuple): (ii,jj)=ii
    a=self.avalue(ii,jj,imap=imap)
    g=self.levs().weight(jj,imap=imap)
    return a*g

# ------------------------------------------------------------------------------

  def gfvalue(self,i,j=None,imap=False):
    """ 
    Return the gf-value for the i-j transition.
    
    Note: it is assumed that state i has a lower energy than state j.
    """
    ii,jj=i,j
    if jj is None and isinstance(ii,tuple): (ii,jj)=ii
    ga=self.gavalue(ii,jj,imap=imap)
    factor=abs(self.levs().energy(ii,imap=imap)-\
               self.levs().energy(jj,imap=imap))
    factor=0.66702*(factor*INVERSECM)**2
    return ga/factor

# ------------------------------------------------------------------------------

  def lifetime(self,i,imap=False):
    """ Return the lifetime of level, i.  Return None if infinite. """
    ii=i
    if imap:
      ii=self.v_levs.lmap(ii)
      if ii is None: raise ErrorRAP('invalid mapping')
    if ii not in self.v_lifetime: return None
    return self.v_lifetime[ii]

# ------------------------------------------------------------------------------

  def prune(self,fraction=0.001):
    """ 
    Remove A-values which compose less than a fraction (default: 0.001) of 
    the (inverse) lifetime of the upper state.
    """
    trans=self.v_data.keys()
    for (i,j) in trans:
      chk=1./self.v_lifetime[j]
      if self.v_data[(i,j)] < fraction*chk: del self.v_data[(i,j)]

# ------------------------------------------------------------------------------

  def set_data(self,i,j=None,value=None,imap=False):
    """ 
    Set data value for the i-j transition.  This is a replacement for the
    set_data() function in rap_trans, as to allow for the lifetime to be 
    properly adjusted.
    
    Note: it is assumed that state i has a lower energy than state j.
    """
    if isinstance(i,tuple):
      (ii,jj)=i
      tval=j
    else:
      ii,jj=i,j
      tval=value
    if tval is None:
      raise ErrorRAP("need to give a new value")
    if ii not in self.levs():
      raise ErrorRAP("Invalid state index given: %s"%ii)
    if jj not in self.levs():
      raise ErrorRAP("Invalid state index given: %s"%jj)
    trans=(ii,jj)
    if imap:
      trans=self.v_levs.lmap(trans)
      if trans is None: raise ErrorRAP('invalid mapping')

    oldval=0.0
    if trans in self: oldval=self[trans]
    self.v_data[trans]=tval

    # need to adjust lifetime
    (ii,jj)=trans
    tau=self.lifetime(jj)
    if tau is None:
      taup=1./tval
    else:
      taup=1./(1./tau-oldval+tval)
    self.v_lifetime[jj]=taup

# ------------------------------------------------------------------------------

  def LS_average(self,omit_incomplete=True):
    """
    Recast LSJ radiative data into LS coupling.  New values are averaged over
    the initial (upper) states and summed over final (lower) states.  The
    option, omit_incomplete, allows the skipping of terms which do not have all
    fine-structure levels present. 
    """
    if not isinstance(self.v_levs,rap_states_lsj):
      raise ErrorRAP("radiative data need to be in LSJ coupling")

    levs_ls=self.v_levs.LS_average(omit_incomplete=omit_incomplete)


    # for each LSJ state, get weighting factors for the constituent levels
    tmap={}        # level -> term map
    wmap={}        # weighting factor for each level
    for idx in levs_ls:
      (conf,term,alfa,par,e)=levs_ls.level(idx)
      wsum=float(rlf_termweight(term))

      tidxs=self.v_levs.search(conf=conf,term=term,alpha=alfa)
      for tidx in tidxs:
        wgt=self.v_levs.weight(tidx)
        wmap[tidx]=float(wgt)/wsum
        tmap[tidx]=idx

    # average rates
    toadd={}
    for (i,j),aval in self.v_data.iteritems():
      t=(tmap[i],tmap[j])
      if t not in toadd: toadd[t]=0.0
      toadd[t]+=wmap[j]*aval

    # new object
    rad=rap_rad(levs=levs_ls)
    rad.v_lifetime=dict([(i,0.0) for i in rad.v_levs.indices()])
    for (i,j),aval in toadd.iteritems():
      rad.v_data[(i,j)]=aval
      rad.v_lifetime[j]+=aval

    for j in rad.v_lifetime.keys():
      asum=rad.v_lifetime[j]
      if asum == 0.0: 
        del rad.v_lifetime[j]
      else:
        rad.v_lifetime[j]=1.0/asum

    return rad

# ------------------------------------------------------------------------------

  def LSJ_split(self,levs_lsj=None):
    """ 
    Split LS energy radiative data in LSJ states.  Need to pass in the LSJ
    level structure (ideally calculated using rap_states_ls.LSJ_split() )
    which will be the level structure associated with the new rap_rad object.
    Dipole allowed transitions will be split according to Eq. 14.50 of Cowan
    (1981).  This assumes that LS coupling is a good description of the system
    and the dipole contribution to the A-value is dominant.  Non-dipole
    transitions will be split statistically and cannot be trusted.  Ideally,
    this routine will serve as a place holder for more accurate values to be
    read in.
    """

    # take term index and get list of level indices (levs_lsj)
    def mapme(i):
      conf=self.levs().configuration(i)
      term=self.levs().term(i)
      alpha=self.levs().seniority(i)
      tlist=levs_lsj.cfilter(conf)
      tlist=levs_lsj.tfilter(term,alpha=alpha,subset=tlist)
      return tlist

    if not isinstance(self.v_levs,rap_states_ls):
      raise ErrorRAP("radiative data need to be in LS coupling")
    if not isinstance(levs_lsj,rap_states_lsj):
      raise ErrorRAP("need to supply rap_states_lsj instance")

    toadd={}
    lmap={}
    for (ilow,iupp) in self:
      if ilow not in lmap: lmap[ilow]=mapme(ilow)
      if iupp not in lmap: lmap[iupp]=mapme(iupp)

      # check if dipole allowed
      dip=True
      (glow,llow)=rlf_termdecode(self.levs().term(ilow))
      (gupp,lupp)=rlf_termdecode(self.levs().term(iupp))
      if glow != gupp: dip=False
      if (llow-lupp) not in [-1,0,+1]: dip=False
      if llow == 0 and lupp == 0: dip=False

      slow=0.5*float(glow-1)
      supp=0.5*float(gupp-1)

      wgts={}
      for ifs_low in lmap[ilow]:
        jlow=levs_lsj.j(ifs_low)

        for ifs_upp in lmap[iupp]:
          jupp=levs_lsj.j(ifs_upp)

          if dip:     # Cowan 14.50
            wgt=(2.*jlow+1.)*(2.*jupp+1.)*\
                pyang_wigner6j(llow,slow,jlow,jupp,1.,lupp)**2
          else:
            wgt=(2.*jlow+1)
          wgts[(ifs_low,ifs_upp)]=wgt

      twgt=sum(wgts.values())
      factor=1./twgt
      aval=self.avalue(ilow,iupp)
      for key in wgts: toadd[key]=factor*wgts[key]*aval

    # assemble new object
    out=rap_rad()
    out.set_levs(levs_lsj)
    out.v_lifetime=dict([(i,0.0) for i in levs_lsj.indices()])
    for (ilow,iupp),aval in toadd.iteritems(): 
      out.v_data[(ilow,iupp)]=aval
      out.v_lifetime[iupp]+=aval
    for j in out.v_lifetime.keys():
      asum=out.v_lifetime[j]
      if asum == 0.0: 
        del out.v_lifetime[j]
      else:
        out.v_lifetime[j]=1.0/asum
    return out

# ------------------------------------------------------------------------------

  def read_olg(self,filename,cup=None,iminlev=None,imaxlev=None,
                                      fminlev=None,fmaxlev=None):
    """ Read radiative data from olg file """
    if self.levs() == None: self.v_levs=rlf_loadstates(filename,cup=cup)
    tcup=self.v_levs.coupling()
    if tcup == 'LSJ':
      count=2
      j0=9
      j1=18
      i0=18
      i1=22
      ae0=33
      ae1=42
      am0=48
      am1=57
    elif tcup == 'LS':
      count=1
      j0=5
      j1=13
      i0=13
      i1=17
      ae0=23
      ae1=32
    else:
      msg='coupling scheme must be LS or LSJ for olg files'
      raise ErrorRAP(msg)

    self.clear(keeplevs=True)
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])
    processing=False
    fin=open(filename,'r')
    for line in fin:
      if 'E0-DATA' in line: count-=1
      if count < 0: break
      if count > 0: continue

      if not processing:
        if '-DATA' in line:
          processing=True
          tmp=re.findall(r'E(\d+).*DATA',line)
          rank=int(tmp[0])
        continue

      if line.strip() == '':
        processing=False
        continue
      j=int(line[j0:j1])
      i=int(line[i0:i1])
      if fminlev and j < fminlev: continue
      if fmaxlev and j > fmaxlev: 
        processing=False
        continue
      try:
        a=float(line[ae0:ae1])     # fails when < 1E-99 (e.g. 6.2456-101)
      except:
        a=0.0
      if tcup == 'LSJ' and rank > 1: a+=float(line[am0:am1])
      self.v_lifetime[j]+=a
      if iminlev and i < iminlev: continue
      if imaxlev and i > imaxlev: continue
      trans=(i,j)
      if trans not in self.v_data: self.v_data[trans]=0.0
      self.v_data[trans]+=a
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def read_adf04(self,filename,cup=None,iminlev=None,imaxlev=None,
                                        fminlev=None,fmaxlev=None):
    """ Read A-values from adf04 file """
    if self.levs() == None: self.v_levs=rlf_loadstates(filename,cup=cup)
    self.clear(keeplevs=True)
    tcup=self.v_levs.coupling()
    if cup is not None and cup != tcup:
      msg='file coupling scheme, %s, incompatible with given: %s'%(tcup,cup)
      raise ErrorRAP(msg)
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])
    fin=open(filename,'r')
    while 1:
      line=fin.readline()
      if line.split()[0] == '-1': break
    line=fin.readline()
    line=fin.readline()
    avalues=[]
    tmp=line.split()
    j=tmp[0]
    i=tmp[1]
    a=tmp[2]
    icol1=line.index(j)+len(j)
    icol2=line[icol1:].index(i)+len(i)+icol1
    icol3=line[icol2:].index(a)+len(a)+icol2
    j=int(j)
    i=int(i)
    a=rlf_adas2float(a)
    if fminlev and j >= fminlev and fmaxlev and j <= fmaxlev:
      self.v_lifetime[j]+=a
      if iminlev and i >= iminlev and imaxlev and i <= imaxlev:
        trans=(i,j)
        self.v_data[trans]=a
    while 1:
      line=fin.readline()
      if line.strip() == '-1': break
      try:
        j=int(line[:icol1])
        i=int(line[icol1:icol2])
      except:
        break
      if (fminlev and j < fminlev) or (fmaxlev and j > fmaxlev): continue
      a=rlf_adas2float(line[icol2:icol3])
      self.v_lifetime[j]+=a
      if (iminlev and i < iminlev) or (imaxlev and i > imaxlev): continue
      trans=(i,j)
      self.v_data[trans]=a
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def read_chianti(self,filename,iminlev=None,imaxlev=None,
                                 fminlev=None,fmaxlev=None):
    """ Read A-values from CHIANTI wfga file """
    id_idx=0        # level index of lower state
    id_jdx=1        # level index of upper state
    id_wl=2         # wavelength (Angstroms)
    id_gf=3         # gf-value (1/s); can be zero if not given
    id_aval=4       # A-value (1/s)
    id_iconf=5      # lower state configuration
    id_istate=6     # lower state term and J-value in spectroscopic notation
    id_dash=7       # separator between upper and lower state information
    id_jconf=8      # upper state configuration
    id_jstate=9     # upper state term and J-value in spectroscopic notation

    if self.levs() == None: raise ErrorRAP('must provide levels object')
    self.clear(keeplevs=True)
    ne=self.levs().num_electrons()
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])

    indices=self.levs().indices()
    fin=open(filename,'r')
    while 1:
      line=fin.readline()
      tmp=line.split()
      if line.strip() == '' or int(tmp[0]) == -1: break
      if len(tmp) != 10: continue
      idx=int(tmp[id_idx])
      jdx=int(tmp[id_jdx])
      if imaxlev and idx > imaxlev: continue
      if fminlev and jdx < fminlev: continue
      if fmaxlev and jdx > fmaxlev: continue
      aval=float(tmp[id_aval])
      iconf=tmp[id_iconf]
      istate=tmp[id_istate]
      jconf=tmp[id_jconf]
      jstate=tmp[id_jstate]
      if idx not in indices or jdx not in indices:
        raise ErrorRAP('level indices not in level structure')

      cchk=self.levs().configuration(idx)
      tchk=self.levs().term(idx)+'%.1f'%(float(0.5*self.levs().j2(idx)))
      ic=rap_config(iconf,ne=ne)
      if ic != str(cchk) or istate != tchk:
        raise ErrorRAP('inconsistency with level labels')

      cchk=self.levs().configuration(jdx)
      tchk=self.levs().term(jdx)+'%.1f'%(float(0.5*self.levs().j2(jdx)))
      jc=rap_config(jconf,ne=ne)
      if jc != str(cchk) or jstate != tchk:
        raise ErrorRAP('inconsistency with level labels')

      self.v_lifetime[jdx]+=aval
      if iminlev and idx < iminlev: continue
      self.v_data[(idx,jdx)]=aval
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def read_tiptop(self,filename,iminlev=None,imaxlev=None,
                                fminlev=None,fmaxlev=None):
    """ 
    Read A-values from TIPTOP gf-value file.

    Note: It is assumed that the symmetry ordering in the TIPTOP files are
    ordered by energy.  That is, the iLV index increases with energy for each
    iSLP symmetry.  If this is not the case, then this routine needs to also
    take in the level file where the symmetries are explicitly defined and
    both the rap_states and rap_rad instances are created here.
    """
    fin=open(filename,'r')

    # get column info
    line=fin.readline()
    hline=fin.readline()
    line=fin.readline()
    tmp=hline.split()
    id_nz=tmp.index('NZ')
    id_ne=tmp.index('NE')
    id_islp=tmp.index('iSLP')
    id_jslp=tmp.index('jSLP')
    id_ilv=tmp.index('iLV')
    id_jlv=tmp.index('jLV')
    id_gf=tmp.index('gF')
    id_wl=tmp.index('WL(A)')

    if self.levs() == None: raise ErrorRAP('must provide levels object')
    self.clear(keeplevs=True)
    ne=self.levs().num_electrons()
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])

    # get mapping from symmetry to level index
    allterms=self.v_levs.terms()
    symmap={}
    for t in allterms:
      (g,l)=rlf_termdecode(t)

      # even
      sym="%s%s0"%(g,l)
      idxs=self.v_levs.symfilter(g,l,0)
      if len(idxs) > 0: symmap[sym]=idxs

      # odd
      sym="%s%s1"%(g,l)
      idxs=self.v_levs.symfilter(g,l,1)
      if len(idxs) > 0: symmap[sym]=idxs

    nz=None
    indices=self.v_levs.indices()
    while 1:
      line=fin.readline()
      tmp=line.split()
      if '=====' in line or line == '': break

      if nz is None:          # check that ion is compatible with levs
        nz=int(tmp[id_nz])
        ne=int(tmp[id_ne])
        if nz != self.v_levs.nuclear_charge() or \
           ne != self.v_levs.num_electrons():
          raise ErrorRAP("incompatible ion")

      islp=tmp[id_islp]
      ilv=int(tmp[id_ilv])
      jslp=tmp[id_jslp]
      jlv=int(tmp[id_jlv])
      gf=float(tmp[id_gf])

      # get indices, want idx to be lower level
      idx=symmap[islp][ilv-1]
      jdx=symmap[jslp][jlv-1]
      if idx not in indices or jdx not in indices:
        raise ErrorRAP('level indices not in level structure')
      if self.v_levs.energy(idx) > self.v_levs.energy(jdx): idx,jdx=jdx,idx
      if imaxlev and idx > imaxlev: continue
      if fminlev and jdx < fminlev: continue
      if fmaxlev and jdx > fmaxlev: continue

      g=self.v_levs.weight(jdx)
      s2=INVERSECM*(self.v_levs.energy(jdx)-self.v_levs.energy(idx))
      s2=s2**2
      aval=GF2A*s2*gf/g

      self.v_lifetime[jdx]+=aval
      if iminlev and idx < iminlev: continue
      self.v_data[(idx,jdx)]=aval
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def read_nist(self,filename,iminlev=None,imaxlev=None,
                              fminlev=None,fmaxlev=None):
    """ Read A-values from NIST text file. """
    id_obwl=0
    id_thwl=1
    id_aval=3
    id_ens=5
    id_confs=6
    id_terms=7
    id_js=8

    if self.levs() == None: raise ErrorRAP('must provide levels object')
    self.clear(keeplevs=True)
    ne=self.levs().num_electrons()
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])

    repchars=['+x','[',']','?']     # these can appear in NIST energy values

    fin=open(filename,'r')
    for line in fin:
      if '---' in line: continue
      if '?' in line: line=line.replace('?',' ')   # these can be annoying
      tmp=line.split('|')
      try:
        obwl=float(tmp[id_obwl])
      except:
        obwl=-1.0
      try:
        thwl=float(tmp[id_thwl])
      except:
        thwl=-1.0
      if obwl < 0.0 and thwl < 0.0: continue
      try:
        aval=float(tmp[id_aval])
      except:
        continue
      if tmp[id_ens].strip() == '': continue
      ttmp=tmp[id_ens].split('-')
      elow=ttmp[0]
      eupp=ttmp[1]
      for ch in repchars:
        elow=elow.replace(ch,'')
        eupp=eupp.replace(ch,'')
      elow=float(elow)/INVERSECM
      eupp=float(eupp)/INVERSECM
      ttmp=tmp[id_confs].split('-')
      clow=rap_config(ttmp[0],ne=ne)
      cupp=rap_config(ttmp[1],ne=ne)
      ttmp=tmp[id_terms].split('-')
      tlow=ttmp[0].replace('*','').strip()
      if '[' in tlow: tlow=None
      tupp=ttmp[1].replace('*','').strip()
      if '[' in tupp: tupp=None
      ttmp=tmp[id_js].split('-')
      tj2low=ttmp[0].strip()
      tj2upp=ttmp[1].strip()
      j2low=int(2*eval(tj2low+'.'))
      j2upp=int(2*eval(tj2upp+'.'))
      idx=self.levs().search(clow,tlow,j2low,energy=elow)
      jdx=self.levs().search(cupp,tupp,j2upp,energy=eupp)

      if fminlev and jdx < fminlev: continue
      if fmaxlev and jdx > fmaxlev: continue
      if imaxlev and idx > imaxlev: continue
      self.v_lifetime[jdx]+=aval
      if iminlev and idx < iminlev: continue
      self.v_data[(idx,jdx)]=aval
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def read_xstar(self,filename,iminlev=None,imaxlev=None,
                               fminlev=None,fmaxlev=None):
    """ Read A-values from xstar type 50 file """
    id_idx=9        # level index of lower state
    id_jdx=10       # level index of upper state
    id_wl=6         # wavelength (Angstroms)
    id_gf=7         # gf-value (1/s); can be zero if not given
    id_aval=8       # A-value (1/s)
    id_ionno=12     # xstar ion number

    if self.levs() == None: raise ErrorRAP('must provide levels object')
    self.clear(keeplevs=True)
    ne=self.levs().num_electrons()
    self.set_initial_limit(iminlev,imaxlev)
    self.set_final_limit(fminlev,fmaxlev)
    self.v_lifetime=dict([(i,0.0) for i in self.v_levs.indices()])

    fin=open(filename,'r')
    while 1:
      line=fin.readline()
      tmp=line.split()
      if line.strip() == '' or int(tmp[0]) == -1: break
      if tmp[-1] != '%' or int(tmp[0]) != 50:
        print line
        raise ErrorRAP('do not understand xstar type 50 file')
      idx=int(tmp[id_idx])
      jdx=int(tmp[id_jdx])
      if idx not in self.v_levs: continue
      if jdx not in self.v_levs: continue
      if fminlev and jdx < fminlev: continue
      if fmaxlev and jdx > fmaxlev: continue
      if imaxlev and idx > imaxlev: continue
      aval=float(tmp[id_aval])
      ionno=int(tmp[id_ionno])
      (z,ne)=xstar_iondecode(ionno)
      if z != self.levs().nuclear_charge() or ne != self.levs().num_electrons():
        raise ErrorRAP('ion in type 50 file does not match level structure')

      self.v_lifetime[jdx]+=aval
      if iminlev and idx > iminlev: continue
      self.v_data[(idx,jdx)]=aval
    fin.close()

    for j in self.v_lifetime.keys():
      asum=self.v_lifetime[j]
      if asum == 0.0: 
        del self.v_lifetime[j]
      else:
        self.v_lifetime[j]=1.0/asum

# ------------------------------------------------------------------------------

  def write_xstar(self,filename,subset=None,precision=5,imap=True):
    """ 
    Write xstar type 50 file. Must have ionization potential and nuclear charge
    specified in levels structure.
    """
    if self.levs().ionization_threshold() == None:
      raise ErrorRAP('must define ionization threshold')
    if self.levs().nuclear_charge() == None:
      raise ErrorRAP('must define nuclear charge')
    indices=self.levs().indices()
    if subset is not None: indices=subset
    trans=[]
    for (i,j) in self.transitions():
      if i in indices and j in indices: trans.append((i,j))

    zeff=self.levs().nuclear_charge()
    ne=self.levs().num_electrons()
    ionno=xstar_ionencode(zeff,ne)

    aform='%'+str(8+precision)+'.'+str(precision)+'E'
    iform='%6i'
    prefix=iform%50+iform%4+iform%0+iform%3+iform%4+iform%0
    suffix=iform%zeff+iform%ionno+' %'
    fout=open(filename,'w')
    for (i,j) in trans:
      try:
        wl=self.wavelength(i,j)
        gf=self.gfvalue(i,j)
      except:
        wl=0.0
        gf=0.0
      av=self.avalue(i,j)
      if av == 0.0: continue
      ii=i
      jj=j
      if imap:
        ii=self.v_levs.lmap_reverse(i)
        jj=self.v_levs.lmap_reverse(j)
        if ii is None or jj is None: continue
      line=prefix+aform%wl+aform%gf+aform%av+iform%ii+iform%jj+suffix
      fout.write(line+'\n')
    fout.close()


################################################################################
class rap_states_jp(rap_states):
  """ Store set of Jp states. """
  def __init__(self,filename=None,**args):
    self.clear()
    if isinstance(filename,rap_states_lsj):
      self._take_lsj(filename)
    elif isinstance(filename,rap_states_jj):
      self._take_jj(filename)
    elif filename: self.loadfile(filename,**args)

# ------------------------------------------------------------------------------

  def _reprocess(self):
    """ Redefine ordinal numbers based on energy-order. """
    if len(self) == 0: return
    indices=zip(*sorted(zip(self.v_energy.values(),self.v_energy.keys())))[1]

    have={}
    tmp={}
    for idx in indices:
      (j2,par,ordi,en)=self.level(idx)

      # get and store seniority
      key=(j2,par)
      if key not in have: have[key]=0
      have[key]+=1
      self.v_ordinal[idx]=have[key]

# ------------------------------------------------------------------------------

  def _take_lsj(self,stlsj):
    """ Take in levels from rap_states_lsj instance. """
    if not isinstance(stlsj,rap_states_lsj):
      msg='expecting rap_states_lsj instance'
      raise ErrorRAP(msg)

    self.clear()
    self.set_ionization_threshold(stlsj.ionization_threshold())
    self.set_core_energy(stlsj.core_energy())
    self.set_nuclear_charge(stlsj.nuclear_charge())
    self.set_num_electrons(stlsj.num_electrons())
    toadd=[]
    for i in stlsj:
      (conf,term,alfa,j2,par,en)=stlsj.level(i)
      toadd.append((en,i,j2,par))
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def _take_jj(self,stjj):
    """ Take in levels from rap_states_jj instance. """
    if not isinstance(stjj,rap_states_jj):
      msg='expecting rap_states_jj instance'
      raise ErrorRAP(msg)

    self.clear()
    self.set_ionization_threshold(stjj.ionization_threshold())
    self.set_core_energy(stjj.core_energy())
    self.set_nuclear_charge(stjj.nuclear_charge())
    self.set_num_electrons(stjj.num_electrons())
    toadd=[]
    for i in stjj:
      (jconf,j2,alfa,par,en)=stjj.level(i)
      toadd.append((en,i,j2,par))
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def _parse_nist_line(self,line):
    """ 
    Function to extract level information (if any) from a NIST data file 
    (from webpage).  
    """
    if 'LIMIT' in line.upper():
      tmp=re.findall(r'(\d+\.\d+)',line)
      if not tmp: return False
      sym=line.split()[0].capitalize()
      zeff=None
      if sym in ATOMSYMBOLS: zeff=ATOMSYMBOLS.index(sym)
      ionth=float(tmp[-1])
      return (zeff,ionth)
    parts=line.split('|')
    if parts[0] == line: return False
    parity=None
    if parts[1].strip() != '':
      parity=0
      if '*' in parts[1]: parity=1
    tmp=re.findall(r'(\d+)/2',parts[2])
    if tmp == []:
      tmp=re.findall(r'(\d+)',parts[2])
      if tmp == []: return False
      jay=float(tmp[0])
    else:
      jay=0.5*float(tmp[0])
    tmp=re.findall(r'(\d+\.\d+)',parts[3])
    energy=None
    if tmp != []: energy=float(tmp[0].replace('+x',''))
    ne=None
    if energy == 0.0:
      conf=rap_config(parts[0])
      ne=conf.num_electrons()
    j2=int(2.0*jay)
    return (j2,parity,energy,ne)

# ------------------------------------------------------------------------------

  def clear(self):
    """ Erase all data. """
    rap_states.clear(self)
    self.v_cup='JP'
    self.v_j2={}
    self.v_ordinal={}

# ------------------------------------------------------------------------------

  def clear_levels(self):
    """ Just erase level data. """
    rap_states.clear_levels(self)
    self.v_cup='JP'
    self.v_j2={}
    self.v_ordinal={}

# ------------------------------------------------------------------------------

  def loadfile(self,filename,**args):
    """ Load in level information from level data file. """
    ft=rlf_filetype(filename)
    if ft == 'adf04':
      self.load_adf04(filename,**args)
    elif ft == 'olg':
      self.load_olg(filename,**args)
    elif ft == 'nist':
      self.load_nist(filename,**args)
    elif ft == 'chianti_elvlc':
      self.load_chianti(filename,**args)
    elif ft == 'x006':
      self.load_xstar(filename,**args)
    elif ft == 'fac_lev':
      self.load_fac(filename,**args)
    else:
      msg='do not know how to read file: %s'%filename
      raise ErrorRAP(msg)

# ------------------------------------------------------------------------------

  def reindex(self,remap=None):
    """ 
    Reassign level indices.  Indices are remapped according to dictionary,
    remap, old (key) -> new (value).  If remap is None, then reindex according
    to energy order.
    """
    imap=rap_states.reindex(self,remap)
    self.v_j2=dict([(imap[i],val) for (i,val) in self.v_j2.items()])
    self.v_parity=dict([(imap[i],val) for (i,val) in self.v_parity.items()])
    self.v_ordinal=dict([(imap[i],val) for (i,val) in self.v_ordinal.items()])

# ------------------------------------------------------------------------------

  def j2range(self,levlist=None,imap=False):
    """ Return tuple containing the minimum and maximum 2J values. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self.mapme,levlist)
    j2s=[self.v_j2[i] for i in levlist if i is not None]
    return (min(j2s),max(j2s))

# ------------------------------------------------------------------------------

  def j2(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_j2[tidx]

# ------------------------------------------------------------------------------

  def jrange(self,levlist=None,imap=False):
    """ Return tuple containing the minimum and maximum 2J values. """
    (j2min,j2max)=self.j2range(levlist,imap)
    return (0.5*float(j2min),0.5*float(j2max))

# ------------------------------------------------------------------------------

  def j(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return 0.5*float(self.v_j2[tidx])

# ------------------------------------------------------------------------------

  def ordinal(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_ordinal[tidx]

# ------------------------------------------------------------------------------

  def level(self,idx,imap=False):
    """ Return all information for single level (j2,par,ordi,E) """
    tidx=self._mapme(idx,imap)
    if tidx is None: return (None,None,None,None)
    out=(self.v_j2[tidx],self.v_parity[tidx],self.v_ordinal[tidx],\
         self.v_energy[tidx])
    return out

# ------------------------------------------------------------------------------

  def label(self,idx,imap=False,fraction=True,spacer=None):
    """ 
    Return level label.

    When fraction is True, return fractional form of J, else decimal form.
    When spacer is not None, use it as an alternative to a space in the label.
    """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    j2=self.j2(idx,imap=imap)
    pval=self.parity(idx,imap=imap)
    ordi=self.ordinal(idx,imap=imap)
    oval=rlf_ordinal(ordi)
    if fraction:
      if j2%2 == 0:
        jstr=str(j2/2)
      else:
        jstr=str(j2)+'/2'
    else:
      jstr='%.1f'%(0.5*float(j2))
    pstr='e'
    if pval == 1: pstr='o'
    out='%s J=%s%s'%(oval,jstr,pstr)
    if spacer is not None: out=out.replace(' ',spacer)
    return out

# ------------------------------------------------------------------------------

  def weight(self,idx,imap=False):
    """ Return weight of level. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.j2(idx,imap=imap)+1

# ------------------------------------------------------------------------------

  def add_level(self,j2,par,energy,idx=False):
    """ Add new level to list """
    if idx is not False and idx in self:
      msg='index already exists: %s'%idx
      raise ErrorRAP(msg)
    tidx=idx
    if tidx is False: 
      tidx=1
      if len(self) > 0: tidx=max(self.indices())+1
    if par not in [0,1]:
      msg='illegal value of parity'
      raise ErrorRAP(msg)
    # check if 2J consistent with other levels (integer v. half-integer)
    if len(self) > 0:
      ichk=self.indices()[0]
      chk1=self.j2(ichk)%2
      chk2=j2%2
      if chk1 != chk2:
        if chk1 == 0:
          msg='new level must have an integral J-value'
        else:
          msg='new level must have a half-integral J-value'
        raise ErrorRAP(msg)
    self.v_j2[tidx]=j2
    self.v_parity[tidx]=par
    self.v_energy[tidx]=energy
    self.v_ordinal[tidx]=-1
    self._reprocess()

# ------------------------------------------------------------------------------

  def add_levels(self,toadd,idxstart=None):
    """ 
    Add set of levels to structure. 

    The input list, toadd, contains tuples with: energy (Ry), index, 2J-value,
    and parity.  If the indexing is to be automated, then set index in the
    tuples to -1.
    """
    if len(toadd) == 0: return
    toadd.sort()
    idxcur=1
    if len(self) > 0: idxcur=max(self.indices())+1
    if idxstart is not None:
      if idxstart < idxcur: 
        msg='invalid starting index: %s'%idxstart
        raise ErrorRAP(msg)
      idxcur=idxstart
    tmp={}
    for (en,idx,j2,par) in toadd:
      tidx=idx
      if tidx in self or tidx in tmp: 
        msg='duplicate index: %s'%tidx
        raise ErrorRAP(msg)
      if tidx < 0:
        tidx=idxcur
        idxcur+=1
      tmp[tidx]=(en,j2,par)
    for idx,(en,j2,par) in tmp.iteritems():
      self.v_j2[idx]=j2
      self.v_parity[idx]=par
      self.v_energy[idx]=en
      self.v_ordinal[idx]=-1
    self._reprocess()

# ------------------------------------------------------------------------------

  def set_num_electrons(self,nelec):
    """ Set number of electrons. """
    self.v_nelec=int(nelec)

# ------------------------------------------------------------------------------

  def set_lmap(self,other,simple=False):
    """ Map to other rap_states_jp instance. """
    if not isinstance(other,rap_states_jp):
      msg='argument must be of type rap_states_jp'
      raise ErrorRAP(msg)

    if simple:
      tmap={}
      for i in other:
        k=None
        for j in self:
          if self[j] == other[i]:
            k=j
            break
        tmap[i]=k
      self.v_lmap=tmap
      self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items()])
      return

    # group level indices by (2J,par)
    fromlist={}
    for i in self:
      key=(self.j2(i),self.parity(i))
      if key not in fromlist: fromlist[key]=[]
      fromlist[key].append(i)
    tolist={}
    for j in other:
      key=(other.j2(j),other.parity(j))
      if key not in tolist: tolist[key]=[]
      tolist[key].append(j)

    self.v_lmap=self._map_by_energy(other,fromlist,tolist)
    for j in other:
      if j not in self.v_lmap: self.v_lmap[j]=None
    self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items() if v is not None])
    for i in self:
      if i not in self.v_lmap_r: self.v_lmap_r[i]=None

# ------------------------------------------------------------------------------

  def j2filter(self,j2list,subset=None,imap=False):
    """ Filter list of states based on 2J values in j2list. """
    if not isinstance(j2list,list) and not isinstance(j2list,tuple): 
      j2list=[j2list]
    tj2list=map(int,j2list)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.j2(i) in tj2list: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def jfilter(self,jlist,subset=None,imap=False):
    """ Filter list of states based on J values in j2list. """
    if not isinstance(jlist,list) and not isinstance(jlist,tuple): jlist=[jlist]
    tjlist=map(float,jlist)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.j(i) in tjlist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def symfilter(self,j2,p=None,subset=None,imap=False):
    """ Filter list of states based on symmetry. """
    if p is None:
      try:
        (jj2,pp)=j2
      except:
        msg='do not understand symmetry'
        raise ErrorRAP(msg)
    else:
      jj2=j2
      pp=p

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      tj2=self.j2(i)
      tp=self.parity(i)
      if tj2 != jj2 or tp != pp: continue
      out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def load_olg(self,filename):
    """ Load in level information from an AUTOSTRUCTURE olg file """
    fin=open(filename,'r')

    # get nuclear charge
    while 1:
      line=fin.readline()
      if line == '':
        msg='file ended before retrieval of nuclear charge'
        raise ErrorRAP(msg)
      if 'ATOMIC NUMBER' in line: break
    tmp=re.findall(r'ATOMIC NUMBER\s*(\d+)',line)
    zeff=int(tmp[0])
    tmp=re.findall(r'NUMBER OF ELECTRONS\s*(\d+)',line)
    ne=int(tmp[0])

    # get level information
    while 1:
      line=fin.readline()
      if line == '': 
        msg='no level data found'
        raise ErrorRAP(msg)
      if '(EK-E1)/RY' in line: break
    tmp=re.findall(r'(-?\d+\.\d+)\s*$',line)
    core_energy=None
    if tmp: core_energy=float(tmp[0])
    toadd=[]
    while 1:
      line=fin.readline()
      tmp=line.split()
      if tmp == []: break
      idx=int(tmp[0])
      j2=int(tmp[6])
      p=0
      if int(tmp[4]) < 0: p=1
      energy=float(tmp[8])
      toadd.append((energy,idx,j2,p))
    fin.close()
    self.set_nuclear_charge(zeff)
    self.set_num_electrons(ne)
    self.set_core_energy(core_energy)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_adf04(self,filename):
    """ Load in level information from an ADAS adf04 file """
    fin=open(filename,'r')

    # basic info
    line=fin.readline()
    tmp=re.findall(r'([\d.]+)',line)
    z=int(tmp[1])
    q=int(tmp[2])
    ip=float(tmp[3])
    ne=z-q+1

    # read in level data
    isic=True
    toadd=[]
    while 1:
      line=fin.readline()
      if '-1' in line: break
      tmp=re.findall(r'(\d+)(.+)\((.+)\)(.+)\((.+)\)\s+(\d+[\.\d]*)',line)[0]
      idx=int(tmp[0])
      conf=rap_config(tmp[1],ne=ne)
      p=conf.parity()
      g=int(tmp[2])
      l=int(tmp[3])
      j2=int(2.*float(tmp[4]))
      e=float(tmp[5])/INVERSECM
      toadd.append((e,idx,j2,p))

      wgt=j2+1
      if wgt not in  range(abs(2*l-g+1)+1,2*l+g+1,2): isic=False
    fin.close()

    # check for correct coupling scheme
    if not isic:
      msg='adf04 file must be LSJ in order to load into rap_states_jp'
      raise ErrorRAP(msg)

    self.set_nuclear_charge(z)
    self.set_num_electrons(ne)
    self.set_ionization_threshold(ip)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_nist(self,filename):
    """ Load in NIST data file (text file taken from website). """
    zeff=None
    nelec=None
    toadd=[]
    ionth=[]
    parity=None
    fin=open(filename,'r')
    for line in fin:
      tmp=self._parse_nist_line(line)
      if tmp:
        if len(tmp) == 2:
          if tmp[0] is not None: zeff=tmp[0]
          ionth.append(tmp[1])
        else:
          (j2,par,en,ne)=tmp
          if par is not None: parity=par
          if ne is not None: nelec=ne
          toadd.append((en,-1,j2,parity))
    fin.close()

    # remove states with no energy
    for i in range(len(toadd)-1,-1,-1):
      if toadd[i][0] is None: toadd.pop(i)

    if len(ionth) > 0: ip=min(ionth)
    self.set_nuclear_charge(zeff)
    if nelec is not None: self.set_num_electrons(nelec)
    self.set_ionization_threshold(ip)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_chianti(self,filename,priority="experimental",only=False):
    """ 
    Read energy levels from CHIANTI file.  CHIANTI energy level files contain
    both theoretical and experimental energies.  The priority option allows
    you to select which energy takes precedence (if available).  The only flag
    will restrict this function to storing only levels where the priority
    energy is present.
    """
    id_idx=0         # level index
    id_k=1           # configuration index
    id_conf=2        # configuration string
    id_g=3           # 2S+1
    id_lval=4        # L (integer)
    id_l=5           # L (string)
    id_j=6           # J
    id_wgt=7         # weight (2J+1)
    id_obcm=8        # observed energy (1/cm)
    id_obry=9        # observed energy (Ry); bad conversion, use 1/cm value
    id_thcm=10       # theoretical energy (1/cm)
    id_thry=11       # theoretical energy (Ry); bad conversion, use 1/cm value

    # attempt to get nuclear charge from filename
    tmp=re.findall(r'(\w+)_\d+\.elvlc',filename)
    if tmp:  
      zeff=ATOMSYMBOLS.index(tmp[0].capitalize())
      self.set_nuclear_charge(zeff)

    ne=None
    toadd=[]
    fin=open(filename,'r')
    for line in fin:
      tmp=line.split()
      if tmp[0] == '-1': break
      if len(tmp) != 12:
        raise RuntimeError('cannot undertand CHIANTI energy level file')
      idx=int(tmp[id_idx])
      conf=rap_config(tmp[id_conf],ne=ne)
      parity=conf.parity()
      g=int(tmp[id_g])
      l=int(tmp[id_lval])
      term=rlf_termencode(g,l)
      j2=int(tmp[id_wgt])-1
      ob=float(tmp[id_obcm])/INVERSECM
      th=float(tmp[id_thcm])/INVERSECM
      if priority[0] == 't':
        en=th
        if en == 0.0 and idx > 1: 
          if only is True: continue
          en=ob
      else:
        en=ob
        if en == 0.0 and idx > 1: 
          if only is True: continue
          en=th
      if ne is None: ne=conf.num_electrons()
      toadd.append((en,idx,j2,parity))
    fin.close()
    self.set_num_electrons(ne)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_xstar(self,filename):
    """ 
    Load level information from xstar type 006 file.
    """
    # column id's
    id_energy=7
    id_j2=8
    id_ip=10
    id_g=12
    id_l=13
    id_z=14
    id_ion=16
    id_id=15
    id_conf=17

    # read level list from file
    ion=None
    ne=None
    ldata=[]
    have=[]
    xid0=None
    fin=open(filename,'r')
    for line in fin:
      data=line.split()

      confstr=data[id_conf-1]
      xid=int(data[id_id-1])
      ion=int(data[id_ion-1])
      g=abs(int(data[id_g-1]))
      l=int(data[id_l-1])
      j2=int(float(data[id_j2-1]))-1
      w=g*(2*l+1)
      if g == 1 or l == 0:
        cup='XX'
      elif w < j2+1:
        cup='CA'
        if 'CA' not in have: have.append('CA')
      elif w == j2+1:
        cup='LS'
        if 'LS' not in have: have.append('LS')
      else:
        j2min=abs(g-1-2*l)
        j2max=g-1+2*l
        okay=False
        for j2t in range(j2min,j2max+1,2):
          if j2t == j2: okay=True
        if okay:
          cup='IC'
          if 'IC' not in have: have.append('IC')
        else:
          cup='??'
      e=float(data[id_energy-1])
      if e == 0.0:
        ne=rap_config(confstr).num_electrons()
        zeff=int(data[id_z-1])
        ion_thresh=float(data[id_ip-1])
      chk=confstr.lower()
      if 'su' in chk or 'sp' in chk: continue 
      if 'co' in chk or 'cn' in chk: continue
      ldata.append([e/RYDBERG,cup,confstr,g,l,j2,xid])
    fin.close()
    if ne is None: raise RuntimeError('cannot find ground state')

    # id's for ldata
    id_e=0
    id_cup=1
    id_conf=2
    id_g=3
    id_l=4
    id_j2=5
    id_xid=6

    # get number of electrons from ground state and standardize configurations
    for i in range(len(ldata)):
      ldata[i][id_conf]=rap_config(ldata[i][id_conf],ne=ne)

    # get total weight of each configuration
    weight={}
    wtot={}
    cups={}
    for i in range(len(ldata)):
      conf=str(ldata[i][id_conf])
      if conf not in wtot: wtot[conf]=ldata[i][id_conf].weight()
      if conf not in cups: cups[conf]=[]
      if ldata[i][id_cup] not in cups[conf]: cups[conf].append(ldata[i][id_cup])
      if ldata[i][id_cup] in ['XX','LS']:
        if conf not in weight:
          weight[conf]=ldata[i][id_g]*(2*ldata[i][id_l]+1)
        else:
          weight[conf]+=ldata[i][id_g]*(2*ldata[i][id_l]+1)
      else:
        if conf not in weight:
          weight[conf]=(ldata[i][id_j2]+1)
        else:
          weight[conf]+=(ldata[i][id_j2]+1)

    # check if configurations have a common coupling scheme
    for c in cups.keys():
      tmp=cups[c][:]
      if 'XX' in tmp: tmp.remove('XX')
      if len(tmp) > 1:
        raise RuntimeError('mixed coupling scheme not supported')
      if tmp == [] and len(have) == 1: tmp=have[:]
      cups[c]=tmp[:]

    # check that 'IC' configurations have all possible levels
    for c in weight.keys():
      if 'IC' not in cups[c]: continue
      elif weight[c] > wtot[c]:
        raise RuntimeError('configuration, '+str(c)+', has too many levels')

    # add levels having IC coupling
    toadd=[]
    for lev in ldata:
      idx=lev[id_xid]
      conf=lev[id_conf]
      sc=str(conf)
      if 'IC' not in cups[sc] and cups[sc] != ['XX']: continue
      j2=lev[id_j2]
      parity=conf.parity()
      energy=lev[id_e]
      toadd.append((energy,idx,j2,parity))
    if len(toadd) == 0:
      msg='no LSJ data found in file, %s'%filename
      raise ErrorRAP(msg)
    self.set_nuclear_charge(zeff)
    self.set_num_electrons(ne)
    self.set_ionization_threshold(ion_thresh/RYDBERG)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_fac(self,filename,nelec=None):
    """ 
    Load in level information from a FAC energy level file. 

    The optional argument, nelec, will tell the function to look for levels with
    the given number of electrons.  Otherwise, the first level set is taken.
    """
    fin=open(filename,'r')

    # basic info
    z=None
    ne=None
    core=None
    nlev=None
    while 1:
      line=fin.readline()
      if line.strip() == '': break
      if ' Z = ' in line:
        z=int(line.split()[-1])
      elif 'E0' in line:
        e0=float(line.split()[-1])

    # loop over blocks
    done=False
    toadd=[]
    while 1:
      if done: break

      # level set info
      while 1:
        line=fin.readline()
        if line == '':
          done=True
          break
        if 'VNL' in line: break
        if 'NELE' in line:
          ne=int(line.split()[-1])
        elif 'NLEV' in line:
          nlev=int(line.split()[-1])
      if nelec is not None and ne != nelec: continue

      # read in level data
      toadd=[]
      while 1:
        line=fin.readline()
        if line.strip() == '': break
        tmp=line.split()
        idx=int(tmp[0])
        en=float(tmp[2])/RYDBERG
        j2=int(tmp[5])

        tmp=re.findall(r'([\w+-]+\(\d+\)\d+)',line)
        tconf=' '.join(tmp)
        conf=rap_jconfig(tconf)
        par=conf.parity()

        toadd.append((en,idx,j2,par))
      if len(toadd) > 0: break
    fin.close()

    if len(toadd) == 0:
      if nelec is None:
        msg="level data found found"
      else:
        pl=''
        if nelec > 1: pl='s'
        msg="level data with %i electron%s not found"%(nelec,pl)
      raise ErrorRAP(msg) 

    self.set_nuclear_charge(z)
    if ne is not None: self.set_num_electrons(ne)
    self.set_core_energy(core)
    self.add_levels(toadd)

################################################################################
class rap_states_jj(rap_states):
  """ Store set of JJ states. """
  def __init__(self,filename=None,**args):
    self.clear()
    if filename: self.loadfile(filename,**args)

# ------------------------------------------------------------------------------

  def _reprocess(self):
    """ Redefine senority-like tags. """
    if len(self) == 0: return
    indices=zip(*sorted(zip(self.v_energy.values(),self.v_energy.keys())))[1]
    allowedterms={}

    have={}
    tmp={}
    for idx in indices:
      (conf,j2,sen,par,en)=self.level(idx)

      key=(conf,j2)
      if key not in have: have[key]=0
      have[key]+=1
      self.v_seniority[idx]=have[key]

# ------------------------------------------------------------------------------

  def _map_jj(self,other,conf=False,alpha=False,common=True,minimize=False,\
              **args):
    """ 
    Set level map to rap_states_jj instance. 

    Setting alpha to True (default) means that the level labels can be trusted
    (conf and common are assumed True).  Otherwise, if conf is True, then 
    configuration is assumed to be a good quantum number and mapping is done
    over (conf,J) pairs; mapping is done in energy order unless minimize is set
    to True when energy differences are minimized (only really necessary
    when number of levels per symmetry is different).  Finally, if conf and
    alpha are False and common is True, then mapping is done over states 
    belonging to the common set of configurations and according to (J,p)
    symmetries.  If conf, alpha, and common are all False, then mapping is
    done over all states according to (J,p).  Setting minimize=True in this 
    case will attempt to minimize energy differences between states.
    """
    if alpha is True:         # seniority numbers are trusted
      tmap={}
      for i in other:
        k=None
        for j in self:
          if self[j] == other[i]:
            k=j
            break
        tmap[i]=k
      self.v_lmap=tmap
      self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items()])
      return

    fromidxs=self.indices()
    toidxs=other.indices()
    if common is True:
      fromconfs=set(self.configurations())
      toconfs=set(other.configurations())
      cconfs=list(fromconfs.intersection(toconfs))
      fromidxs=self.cfilter(cconfs)
      toidxs=other.cfilter(cconfs)

    if conf is True:
      fromlist={}
      for i in fromidxs:
        key=(self.configuration(i),self.j2(i))
        if key not in fromlist: fromlist[key]=[]
        fromlist[key].append(i)
      tolist={}
      for j in toidxs:
        key=(other.configuration(j),other.j2(j))
        if key not in fromlist: continue
        if key not in tolist: tolist[key]=[]
        tolist[key].append(j)
    else:
      fromlist={}
      for i in fromidxs:
        key=(self.j2(i),self.parity(i))
        if key not in fromlist: fromlist[key]=[]
        fromlist[key].append(i)
      tolist={}
      for j in toidxs:
        key=(other.j2(j),other.parity(j))
        if key not in fromlist: continue
        if key not in tolist: tolist[key]=[]
        tolist[key].append(j)

    if minimize is True:
      self.v_lmap=self._map_by_energy(other,fromlist,tolist)
    else:
      self.v_lmap=self._map_straight(other,fromlist,tolist)
    for j in other:
      if j not in self.v_lmap: self.v_lmap[j]=None
    self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items() if v is not None])
    for i in self:
      if i not in self.v_lmap_r: self.v_lmap_r[i]=None

# ------------------------------------------------------------------------------

  def _map_lsj(self,other,conflist=None,**args):
    """ 
    Set level map to rap_states_lsj instance. 

    Mapping is performed only over the levels sharing a set of common
    configurations.  The mapping is based on J, parity, and energy order.
    """
    if conflist is not None:
      if isinstance(conflist[0],rap_config):
        oconfs=set(conflist[:])
        mejconfs=set(self.configurations())
      elif isinstance(conflist[0],rap_jconfig):
        oconfs=set(other.configurations())
        mejconfs=set(conflist[:])
      else:
        msg='confset must contains instances of rap_config or rap_jconfig'
        raise ErrorRAP(msg)
    else:
      oconfs=set(other.configurations())
      mejconfs=set(self.configurations())

    meconfs=[]
    for jconf in mejconfs:
      conf=jconf.to_config()
      if conf not in meconfs: meconfs.append(conf)
    meconfs=set(meconfs)
    shared=list(oconfs & meconfs)
    jshared=[]
    for conf in shared: jshared.extend(conf.to_jconfig())
    oidxs=other.cfilter(shared)
    meidxs=self.cfilter(jshared)

    (j2min,j2max)=other.j2range(levlist=oidxs)
    pars=[0,1]
    fromlist={}
    tolist={}
    for p in pars:
      for j2 in range(j2min,j2max+1,2):
        sym=(j2,p)
        toidxs=sorted(other.symfilter(sym,subset=oidxs))
        tmeidxs=sorted(self.symfilter(sym,subset=meidxs))
        if len(toidxs) != len(tmeidxs):        # sanity check
          msg='level sets differ in number of levels for symmetry: '+\
              '2J,p=%i,%i'%(j2,p)
          raise ErrorRAP(msg)
        fromlist[sym]=tmeidxs
        tolist[sym]=toidxs

    self.v_lmap=self._map_by_energy(other,fromlist,tolist)
    for j in other:
      if j not in self.v_lmap: self.v_lmap[j]=None
    self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items() if v is not None])
    for i in self:
      if i not in self.v_lmap_r: self.v_lmap_r[i]=None

# ------------------------------------------------------------------------------

  def _map_jp(self,other,conflist=None,**args):
    """ 
    Set level map to rap_states_jp instance. 

    Mapping is performed only over all levels and based on J, parity, and 
    energy order.
    """
    idxs=self.indices()
    if conflist is not None:
      if isinstance(conflist[0],rap_config):
        tconflist=[]
        for conf in conflist: tconflist.extend(conf.to_jconfig())
        conflist=tconflist
      idxs=self.cfilter(conflist)
    (j2min,j2max)=self.j2range(levlist=idxs)
    pars=[0,1]
    fromlist={}
    tolist={}
    for p in pars:
      for j2 in range(j2min,j2max+1,2):
        sym=(j2,p)
        fromlist[sym]=self.symfilter(sym)
        tolist[sym]=other.symfilter(sym)

    self.v_lmap=self._map_by_energy(other,fromlist,tolist)
    for j in other:
      if j not in self.v_lmap: self.v_lmap[j]=None
    self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items() if v is not None])
    for i in self:
      if i not in self.v_lmap_r: self.v_lmap_r[i]=None

# ------------------------------------------------------------------------------

  def clear(self):
    """ Erase all data. """
    rap_states.clear(self)
    self.v_cup='JJ'
    self.v_jconf={}
    self.v_j2={}
    self.v_seniority={}

# ------------------------------------------------------------------------------

  def clear_levels(self):
    """ Just erase level data. """
    rap_states.clear_levels(self)
    self.v_cup='JJ'
    self.v_jconf={}
    self.v_j2={}
    self.v_seniority={}

# ------------------------------------------------------------------------------

  def loadfile(self,filename,**args):
    """ Load in level information from level data file. """
    ft=rlf_filetype(filename)
    if ft == 'fac_lev':
      self.load_fac(filename,**args)
    else:
      msg='do not know how to read file: %s'%filename
      raise ErrorRAP(msg)

# ------------------------------------------------------------------------------

  def reindex(self,remap=None):
    """ 
    Reassign level indices.  Indices are remapped according to dictionary,
    remap, old (key) -> new (value).  If remap is None, then reindex according
    to energy order.
    """
    imap=rap_states.reindex(self,remap)
    self.v_jconf=dict([(imap[i],val) for (i,val) in self.v_jconf.items()])
    self.v_j2=dict([(imap[i],val) for (i,val) in self.v_j2.items()])
    self.v_seniority=dict([(imap[i],val) for (i,val) in \
                          self.v_seniority.items()])

# ------------------------------------------------------------------------------

  def configurations(self,levlist=None,imap=False):
    """ Return list of all configurations. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self._mapme,levlist)
    out=[]
    for i in levlist:
      if self.v_jconf[i] in out: continue
      out.append(self.v_jconf[i])
    return out

# ------------------------------------------------------------------------------

  def configuration(self,idx,imap=False):
    """ Return configuration of level with given index. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_jconf[tidx]

# ------------------------------------------------------------------------------

  def seniority(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_seniority[tidx]

# ------------------------------------------------------------------------------

  def j2range(self,levlist=None,imap=False):
    """ Return tuple containing the minimum and maximum 2J values. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self.mapme,levlist)
    j2s=[self.v_j2[i] for i in levlist if i is not None]
    return (min(j2s),max(j2s))

# ------------------------------------------------------------------------------

  def j2(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_j2[tidx]

# ------------------------------------------------------------------------------

  def jrange(self,levlist=None,imap=False):
    """ Return tuple containing the minimum and maximum 2J values. """
    (j2min,j2max)=self.j2range(levlist,imap)
    return (0.5*float(j2min),0.5*float(j2max))

# ------------------------------------------------------------------------------

  def j(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return 0.5*float(self.v_j2[tidx])

# ------------------------------------------------------------------------------

  def level(self,idx,imap=False):
    """ Return all information for single level (conf,j2,alpha,par,E) """
    tidx=self._mapme(idx,imap)
    if tidx is None: return (None,None,None,None,None)
    out=(self.v_jconf[tidx],self.v_j2[tidx],self.v_seniority[tidx],\
         self.v_parity[tidx], self.v_energy[tidx])
    return out

# ------------------------------------------------------------------------------

  def label(self,idx,imap=False,spacer=None):
    """ Return level label; optional spacer can be used instead of space. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    j2=self.j2(idx,imap=imap)
    if j2%2 == 0:
      jstr=str(j2/2)
    else:
      jstr=str(j2)+'/2'
    out=str(self.v_jconf[tidx])+' J=%s'%jstr
    alfa=self.v_seniority[tidx]
    if alfa > 1: out+='#%s'%alfa
    if spacer is not None: out=out.replace(' ',spacer)
    return out

# ------------------------------------------------------------------------------

  def weight(self,idx,imap=False):
    """ Return weight of level. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.j2(idx,imap=imap)+1

# ------------------------------------------------------------------------------

  def add_level(self,conf,j2,energy,idx=False):
    """ Add new level to list """
    if idx is not False and idx in self:
      msg='index already exists: %s'%idx
      raise ErrorRAP(msg)
    tidx=idx
    if tidx is False: 
      tidx=1
      if len(self) > 0: tidx=max(self.indices())+1
    tconf=conf
    if not isinstance(tconf,rap_jconfig): tconf=rap_jconfig(tconf)
    if self.num_electrons() is None:
      self.v_nelec=tconf.num_electrons()
    else:
      if tconf.num_electrons() != self.num_electrons():
        msg='inconsistent number of electrons; expect/have: %s, %s'%\
            (self.num_electrons(),tconf.num_electrons())
        raise ErrorRAP(msg)
    self.v_jconf[tidx]=tconf
    self.v_j2[tidx]=j2
    self.v_seniority[tidx]=-1
    self.v_parity[tidx]=tconf.parity()
    self.v_energy[tidx]=energy
    self._reprocess()

# ------------------------------------------------------------------------------

  def add_levels(self,toadd,idxstart=None):
    """ 
    Add set of levels to structure. 

    The input list, toadd, contains tuples with: energy (Ry), index,
    configuration, and 2J.  If the indexing is to be automated, then set index 
    in the tuples to -1.
    """
    if len(toadd) == 0: return
    toadd.sort()
    idxcur=1
    if len(self) > 0: idxcur=max(self.indices())+1
    if idxstart is not None:
      if idxstart < idxcur: 
        msg='invalid starting index: %s'%idxstart
        raise ErrorRAP(msg)
      idxcur=idxstart
    tmp={}
    nechk=None
    for (en,idx,conf,j2) in toadd:
      tidx=idx
      if tidx in self or tidx in tmp: 
        msg='duplicate index: %s'%tidx
        raise ErrorRAP(msg)
      if tidx < 0:
        tidx=idxcur
        idxcur+=1
      tconf=conf
      if not isinstance(conf,rap_jconfig): tconf=rap_jconfig(conf)
      if nechk is None: nechk=tconf.num_electrons()
      if tconf.num_electrons() != nechk:
        msg='inconsistent number of electrons; expected/have: %s, %s'%\
            (nechk,tconf.num_electrons())
        raise ErrorRAP(msg)
      tmp[tidx]=(en,tconf,j2)
    self.v_nelec=nechk
    for idx,(en,conf,j2) in tmp.iteritems():
      self.v_jconf[idx]=conf
      self.v_j2[idx]=j2
      self.v_seniority[idx]=-1
      self.v_parity[idx]=conf.parity()
      self.v_energy[idx]=en
    self._reprocess()

# ------------------------------------------------------------------------------

  def set_lmap(self,other,**args):
    """ Map to other rap_states_jj instance. """
    if self.num_electrons() != other.num_electrons():
      msg='inconsistent number of electrons; expected/have: %s, %s'%\
          (self.num_electrons(),other.num_electrons())
      raise ErrorRAP(msg)
    if isinstance(other,rap_states_jj):
      self._map_jj(other,**args)
    elif isinstance(other,rap_states_lsj):
      self._map_lsj(other,**args)
    elif isinstance(other,rap_states_jp):
      self._map_jp(other,**args)
    else:
      msg='argument must be of type rap_states_jj or rap_states_lsj'
      raise ErrorRAP(msg)

# ------------------------------------------------------------------------------

  def cfilter(self,clist,subset=None,imap=False):
    """ Filter list of states based on configurations in clist. """
    if not isinstance(clist,list) and not isinstance(clist,tuple): clist=[clist]
    tclist=[]
    for c in clist:
      tc=c
      if not isinstance(tc,rap_jconfig): tc=rap_jconfig(tc)
      tclist.append(tc)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.configuration(i) in tclist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def j2filter(self,j2list,subset=None,imap=False):
    """ Filter list of states based on 2J values in j2list. """
    if not isinstance(j2list,list) and not isinstance(j2list,tuple): 
      j2list=[j2list]
    tj2list=map(int,j2list)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.j2(i) in tj2list: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def jfilter(self,jlist,subset=None,imap=False):
    """ Filter list of states based on J values in j2list. """
    if not isinstance(jlist,list) and not isinstance(jlist,tuple): jlist=[jlist]
    tjlist=map(float,jlist)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.j(i) in tjlist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def symfilter(self,j2,p=None,subset=None,imap=False):
    """ Filter list of states based on symmetry. """
    if p is None:
      try:
        (jj2,pp)=j2
      except:
        msg='do not understand symmetry'
        raise ErrorRAP(msg)
    else:
      jj2=j2
      pp=p

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      tj2=self.j2(i)
      tp=self.parity(i)
      if tj2 != jj2 or tp != pp: continue
      out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def search(self,conf=None,j2=None,parity=None,alpha=None,energy=None):
    """ 
    Look for matching levels in structure based on configuration and
    2J values.  If duplicate terms are present and alpha is not given, a list
    of all matching indices are returned.  If alpha is given, the index matching
    that seniority is returned.  If energy is given (and alpha is not), then
    the senority closest to that energy is returned.  If conf and term are 
    both None, then return matching index based on j2, parity, and energy order
    given by alpha.
    """
    if j2 is None: 
      if conf is None:
        msg="if no 2J is specified, must provide configuration"
        raise ErrorRAP(msg)
      levs=self.cfilter(conf)
      return levs

    levs=self.j2filter(j2)
    if conf is None or term is None:
      if conf is not None: parity=conf.parity()
      if parity is None: 
        raise RuntimeError('need parity when conf is None')
      levs=self.pfilter(parity,subset=levs)
      if levs == []: raise RuntimeError('matching level not found')
      if alpha is not None:
        return levs[alpha-1]
      elif energy is not None:
        tmp=[]
        for i in levs: tmp.append((abs(energy-self.energy(i)),i))
        return min(tmp)[1]
      return levs

    levs=self.cfilter(conf,subset=levs)
    if levs == []: raise RuntimeError('matching level not found')
    if len(levs) == 1: return levs[0]
    if alpha is not None:
      for i in levs: 
        if self.seniority(i) == alpha: return i
      raise RuntimeError('matching level not found')
    elif energy is not None:
      tmp=[]
      for i in levs: tmp.append((abs(energy-self.energy(i)),i))
      return min(tmp)[1]
    return levs

# ------------------------------------------------------------------------------

  def load_fac(self,filename,nelec=None):
    """ 
    Load in level information from a FAC energy level file. 

    The optional argument, nelec, will tell the function to look for levels with
    the given number of electrons.  Otherwise, the first level set is taken.
    """
    fin=open(filename,'r')

    # basic info
    z=None
    ne=None
    core=None
    nlev=None
    while 1:
      line=fin.readline()
      if line.strip() == '': break
      if ' Z = ' in line:
        z=int(line.split()[-1])
      elif 'E0' in line:
        e0=float(line.split()[-1])

    # loop over blocks
    done=False
    toadd=[]
    while 1:
      if done: break

      # level set info
      while 1:
        line=fin.readline()
        if line == '':
          done=True
          break
        if 'VNL' in line: break
        if 'NELE' in line:
          ne=int(line.split()[-1])
        elif 'NLEV' in line:
          nlev=int(line.split()[-1])
      if nelec is not None and ne != nelec: continue

      # read in level data
      while 1:
        line=fin.readline()
        if line.strip() == '': break
        tmp=line.split()
        idx=int(tmp[0])
        en=float(tmp[2])/RYDBERG
        j2=int(tmp[5])

        # sub-shell occupancy
        nocc={}
        tmp=re.findall(r'(\d+)\*(\d+)',line)
        for t in tmp: nocc[int(t[0])]=int(t[1])
        ssocc={}
        tmp=re.findall(r'(\d+[a-z])(\d+)',line)
        for t in tmp: ssocc[t[0]]=int(t[1])
        for n,wn in nocc.iteritems():
          wdiff=wn
          noell=True
          for l in range(n):
            ss=rlf_ssencode(n,l)
            if ss not in ssocc: continue
            noell=False
            wdiff-=ssocc[ss]
          if noell:              # one or all sub-shells filled
            if 2*n**2 == wn:
              for l in range(n):
                ss=rlf_ssencode(n,l)
                ssocc[ss]=4*l+2
              continue
            for l in range(n):
              ss=rlf_ssencode(n,l)
              tw=rlf_ssoccupancy(ss)
              done=False
              if not done and tw == wn:   # fill first available
                ssocc[ss]=wn
              else:
                ssocc[ss]=0
            continue
          if wdiff > 0:
            for l in range(n):
              ss=rlf_ssencode(n,l)
              if ss in ssocc: continue
              if rlf_ssoccupancy(ss) <= wdiff:
                ssocc[ss]=rlf_ssoccupancy(ss)
                wdiff-=ssocc[ss]
          if wdiff != 0:                    # sanity check
            msg='inconsistancy in FAC level file?  idx:%i'%idx
            raise ErrorRAP(msg)

        tmp=re.findall(r'([\w+-]+\(\d+\)\d+)',line)
        tconf=' '.join(tmp)
        conf=rap_jconfig(tconf,ne=ne,ssocc=ssocc)

        toadd.append((en,idx,conf,j2))
      if len(toadd) > 0: break
    fin.close()

    if len(toadd) == 0:
      if nelec is None:
        msg="level data found found"
      else:
        pl=''
        if nelec > 1: pl='s'
        msg="level data with %i electron%s not found"%(nelec,pl)
      raise ErrorRAP(msg) 

    self.set_nuclear_charge(z)
    self.set_core_energy(core)
    self.add_levels(toadd)

################################################################################

class rap_config:
  """ Stores a single configuration. """
  def __init__(self,conforne,ne=None):
    if isinstance(conforne,rap_jconfig):
      tmp=conforne.to_config()
      self.v_conf=tmp.v_conf
    elif isinstance(conforne,int):
      self.v_conf=self._baseconfig(conforne)
    else:
      self.v_conf=self._format_config(conforne,ne=ne)

# ------------------------------------------------------------------------------

  def __repr__(self):
    return self.label()

# ------------------------------------------------------------------------------

  def __eq__(self,conf):
    tconf=conf
    if not isinstance(tconf,rap_config): tconf=rap_config(tconf)
    if repr(self) == repr(tconf): return True
    return False

# ------------------------------------------------------------------------------

  def __ne__(self,conf):
    tconf=conf
    if not isinstance(tconf,rap_config): tconf=rap_config(tconf)
    if repr(self) != repr(tconf): return True
    return False

# ------------------------------------------------------------------------------

  def __cmp__(self,conf):
    return rlf_confsorter(self,conf)

# ------------------------------------------------------------------------------

  def __hash__(self):
    return hash(self.label())

# ------------------------------------------------------------------------------

  def __len__(self):
    """ Return length of configuration label. """
    return len(self.label())

# ------------------------------------------------------------------------------

  def __add__(self,other):
    return self.label()+other

# ------------------------------------------------------------------------------

  def __radd__(self,other):
    return other+self.label()

# ------------------------------------------------------------------------------

  def __getitem__(self,idx):
    return self.occupancy(idx)

# ------------------------------------------------------------------------------

  def __iter__(self):
    return self.subshells().__iter__()

# ------------------------------------------------------------------------------

  def _baseconfig(self,ne):
    """ Fill sub-shells in order until number of electrons is reached. """
    netot=ne
    out={}
    i=1
    while netot > 0:
      w=rlf_ssoccupancy(i)
      out[i]=min(w,netot)
      netot-=w
      i+=1
    return out

# ------------------------------------------------------------------------------

  def _format_config(self,inconf,ne=None):
    """ Attempt to understand input configuration. """

    # function to see if configuration is in Eissner form
    def read_eissner(x):
      x=str(x)
      tmp=re.findall(r'(\d+)([a-z])',x)
      if tmp: return False
      if len(x)%3 == 1: return False
      if len(x)%3 == 2: x='5'+x
      out={}
      for i in range(len(x)/3):
        xx=x[3*i:3*i+3]
        try:
          w=int(xx[:2])-50
          ss=rlf_eissner2ss(xx[-1],index=False)
        except:
          return False
        if w > 1 and rlf_ssindex(ss) > 9: return False   # probably not Eissner
        if w < 0 or w > rlf_ssoccupancy(ss): return False
        out[ss]=w
      conf=[]
      for ss,w in out.iteritems(): conf.append(ss+str(w))
      conf=' '.join(conf)
      return conf


    # configuration must be a string
    if isinstance(inconf,rap_config):
      return inconf
    if not isinstance(inconf,str): 
      msg='configuration must be a string; conf: %s'%inconf
      raise ErrorRAP(msg)

    conf=inconf+' '
    out={}

    # look for Eissner notation
    tmp=re.split(r'(\w+\.)',conf.strip())
    while tmp[0] == '': tmp.pop(0)
    chk=read_eissner(tmp[0])
    if chk is not False: conf=chk

    # configuration must contain both numbers and letters
    tmp=re.findall(r'(\d)',conf)
    if tmp == []: 
      msg='configuration contains no digits; conf: %s'%conf
      raise ErrorRAP(msg)
    tmp=re.findall(r'([A-Za-z])',conf)
    if tmp == []: 
      msg='configuration contains no letters; conf: %s'%conf
      raise ErrorRAP(msg)

    # if only uppercase letters exist, assume no term labels
    tmp=re.findall(r'([a-z])',conf)
    if tmp == []: conf=conf.lower()

    # get rid of carets
    conf=conf.replace('^','')

    # remove parity characters
    conf=conf.replace('e','')
    conf=conf.replace('o','')

    # identify any term labels and remove
    tmp=re.findall(r'(#\d+)',conf)
    for t in tmp: conf=conf.replace(t,'')
    tmp=re.findall(r'(\(?\s*\d+[A-Z]_?\d+\.5[^a-z]\s*\)?)',conf)
    for t in tmp: conf=conf.replace(t,'')
    tmp=re.findall(r'(\(?\s*\d+[A-Z]_?[\d/]*\s*\)?)',conf)
    for t in tmp: conf=conf.replace(t,'')
    conf=conf.strip()+' '

    # deal with hole-notation
    sshole=None
    whole=None
    idxhole=None
    tmp=re.findall(r'(\[\d+[a-z][\d]*\])',conf)
    if tmp != []:
      if len(tmp) > 1:
        msg='only one hole supported; conf: %s'%conf
        raise ErrorRAP(msg)
      t=tmp[0]
      ttmp=re.findall(r'(\d+)([a-z])(\d*)',t)
      n=int(ttmp[0][0])
      l=SPECLETTERS.index(ttmp[0][1].upper())
      w=1
      if ttmp[0][2].strip() != '': w=int(ttmp[0][2])
      wmax=4*l+2
      sshole=str(n)+SPECLETTERS[l].lower()
      whole=wmax-w
      if whole < 0: 
        msg='too many electrons in hole; conf: %s'%conf
        raise ErrorRAP(msg)
      conf=conf.replace(t,'')
      idxhole=rlf_ssindex(sshole)

    # check for letters and numbers
    tmp=re.findall(r'(\d)',conf)
    if len(tmp) == 0:
      msg='missing numbers in main part of configuration: %s'%inconf
      raise ErrorRAP(msg)
    tmp=re.findall(r'([A-Za-z])',conf)
    if len(tmp) == 0:
      msg='missing letters in main part of configuration: %s'%inconf
      raise ErrorRAP(msg)

    # need to separate sub-shells; first check for division characters
    # (space, period, etc) and if unsuccessful, do it the hard way
    tmp=re.split(r'[\W\s._]+',conf)
    while tmp[-1] == '': tmp.pop()
    while tmp[0] == '': tmp.pop(0)
    hardway=False
    ttmp=re.findall(r'([a-z])',conf)
    if len(ttmp) != len(tmp): hardway=True
    if not hardway:
      sslist=[]
      for t in tmp:
        if t == '': continue
        ttmp=list(re.findall(r'(\d+)([a-z])(\d*)',t)[0])
        if ttmp == []: 
          msg='cannot resolve sub-shells; conf: %s'%conf
          raise ErrorRAP(msg)
        if ttmp[-1] == '': ttmp[-1]=1
        sslist.append(ttmp)      
    else:

      # get rid of all unnecessary characters
      tmp=re.findall(r'([a-z\d]+)',conf)
      conf=''
      for t in tmp: conf+=t

      # make temporary configuration structure
      ss=re.findall(r'([a-z])',conf)
      sslist=[]
      for s in ss: sslist.append(['',s,''])
      nss=len(sslist)

      # separate digits by choosing largest allowed occupation
      digits=re.findall(r'(\d+)',conf)
      sslist[0][0]=int(digits[0])     # first n value is easy
      for i in range(nss-1):
        s=sslist[i][1]
        l=SPECLETTERS.index(s.upper())
        wmax=4*l+2
        j=len(digits[i+1])-1
        while j > 0 and int(digits[i+1][:j]) > wmax: j-=1
        w=digits[i+1][:j]
        if w == '': w=1
        sslist[i][2]=int(w)
        sslist[i+1][0]=int(digits[i+1][j:])
      sslist[-1][2]=1
      if len(digits) == nss+1: sslist[-1][2]=int(digits[-1])

    # assemble sub-shell dictionary
    ssdict={}
    wmax={}
    for ss in sslist: 
      s=str(ss[0])+ss[1]
      w=int(ss[2])
      if s not in ssdict: ssdict[s]=0
      ssdict[s]+=w
      wmax[s]=4*SPECLETTERS.index(ss[1].upper())+2
      if ssdict[s] > wmax[s]: 
        msg='maximum occupancy exceeded; conf: %s'%conf
        raise ErrorRAP(msg)
    sss=ssdict.keys()
    sss.sort(rlf_sssorter)
    first=rlf_ssindex(sss[0])
    if idxhole and idxhole >= first:
      if sshole not in ssdict: ssdict[sshole]=0
      ssdict[sshole]+=whole
      if sshole in wmax and ssdict[sshole] > wmax[sshole]: 
        msg='maximum occupancy exceeded; conf: %s'%conf
        raise ErrorRAP(msg)
      sss=ssdict.keys()
      sss.sort(rlf_sssorter)
    tconf=[]
    for s in sss: tconf.append((s,ssdict[s]))

    # deal with hole-notation
    idxfillin=[]
    if idxhole and idxhole < first:
      ssdict[sshole]=whole
      for i in range(idxhole+1,first):
        ss=rlf_subshell(i)
        ssdict[ss]=rlf_ssoccupancy(ss)
        idxfillin.append(rlf_ssindex(ss))
    sss=ssdict.keys()
    sss.sort(rlf_sssorter)

    # fill in full sub-shells
    first=rlf_ssindex(sss[0])
    for i in range(1,first):
      ss=rlf_subshell(i)
      ssdict[ss]=rlf_ssoccupancy(ss)
      idxfillin.append(rlf_ssindex(ss))

    # prepare output
    out={}
    for ss in ssdict: 
      if ssdict[ss] > 0: out[rlf_ssindex(ss)]=ssdict[ss]

    # check if configuration has the right number of electrons
    if ne is None: return out
    nediff=sum(ssdict.values())-ne
    if nediff == 0: return out
    if nediff < 0:
      msg='specified number of electrons exceeded; conf: %s'%conf
      raise ErrorRAP(msg)
    if nediff%2 == 1: 
      msg='cannot satisfy given number of electrons; conf: %s'%conf
      raise ErrorRAP(msg)

    # look for single sub-shell outside of a core
    tmp=[ssdict[ss] for ss in sss]
    if tmp.count(0) == len(tmp)-1:
      for ss in sss:
        if ssdict[ss] > 0: break
      nechk=ne-ssdict[ss]
      tout={}
      tout[rlf_ssindex(ss)]=ssdict[ss]
      if nechk == 0: return tout
      i=1
      tot=0
      okay=False
      while 1:
        w=rlf_ssoccupancy(i)
        tot+=w
        tout[i]=w
        if tot == nechk: okay=True
        if tot >= nechk: break
        i+=1
      if okay: return tout

    # search for unspecified empty sub-shells using a combinatorial
    # approach.  select matches based on fewest sub-shells and
    # highest set of sub-shells (2s > 1s and 2s 2p > 1s 2p).
    idxfillin.sort()
    idxfillin=tuple(idxfillin)
    wlist={}
    for idx in idxfillin:
      wlist[idx]=rlf_ssoccupancy(idx)
    match=False
    for n in range(1,len(idxfillin)+1):
      for comb in iter_combinations(idxfillin,n):
        wtot=0
        for idx in comb: wtot+=wlist[idx]
        if wtot == nediff: 
          rank=sum(comb)
          if match is False or rank > sum(match): match=comb
      if match is not False: break
    if match is False: 
      msg='cannot satisfy given number of electrons; conf: %s'%conf
      raise ErrorRAP(msg)
    idxmin=min(match)
    idxmax=max(match)
    for i in range(idxmin,idxmax+1):
      if i in match: del out[i]
    return out

# ------------------------------------------------------------------------------

  def copy(self):
    """ Return a copy. """
    return copy.deepcopy(self)

# ------------------------------------------------------------------------------

  def label(self,showclosed=False):
    if self.num_electrons() == 0: return ''
    out=[]
    if showclosed:
      for ss in self.subshells(): out.append(ss+str(self.occupancy(ss)))
    else:
      imax=max(self.v_conf.keys())
      icore=1
      while icore in self.v_conf and \
            self.v_conf[icore] == rlf_ssoccupancy(icore):
        icore+=1
      if icore > imax: 
        icore=imax
      elif icore not in self.v_conf:
        out.append(rlf_subshell(icore)+'0')
      for i in self.subshells(indices=True):
        w=self.occupancy(i)
        if i < icore: continue
        out.append(rlf_subshell(i)+str(w))
    out=' '.join(out)
    return out

# ------------------------------------------------------------------------------

  def num_electrons(self):
    """ Return total number of electrons. """
    return sum(self.v_conf.values())

# ------------------------------------------------------------------------------

  def subshells(self,indices=False):
    """ Return sorted list of occupied sub-shells. """
    ss=self.v_conf.keys()
    ss.sort()
    if indices: return ss
    return map(rlf_subshell,ss)

# ------------------------------------------------------------------------------

  def occupancy(self,ss):
    """ 
    Return occupancy of sub-shell, ss (which can be index or string form. 
    """
    sss=ss
    if not isinstance(sss,int): sss=rlf_ssindex(sss)
    if sss not in self.v_conf: return 0
    return self.v_conf[sss]

# ------------------------------------------------------------------------------

  def parity(self):
    """ Return parity of configuration. """
    p=1
    for i in self.subshells():
      (n,l)=rlf_ssdecode(i)
      if l%2 == 0: continue
      p*=(-1)**self.occupancy(i)
    return (1-p)/2

# ------------------------------------------------------------------------------

  def weight(self):
    """ Return statistical weight of configuration. """
    wgt=1
    for ss in self.subshells():
      w=self.occupancy(ss)
      (n,l)=rlf_ssdecode(ss)
      wgt*=rlf_binom(4*l+2,w)
    return int(wgt)

# ------------------------------------------------------------------------------

  def valence_subshell(self):
    """ Return valence sub-shell. """
    tss=max(self.subshells(indices=True))
    return rlf_subshell(tss)

# ------------------------------------------------------------------------------

  def pop_electron(self,ss=None):
    """ 
    Return a copy of the configuration where an electron from the given 
    sub-shell has been removed (last sub-shell by default).
    """
    out=copy.deepcopy(self)
    tss=ss
    if tss is None: 
      tss=max(out.subshells(indices=True))
    else:
      if not isinstance(tss,int): tss=rlf_ssindex(tss)
    if tss not in out.subshells(indices=True): 
      msg='specified sub-shell not present: %s'%ss
      raise ErrorRAP(msg)
    out.v_conf[tss]-=1
    if out.v_conf[tss] == 0: del out.v_conf[tss]
    return out

# ------------------------------------------------------------------------------

  def add_electron(self,ss=None):
    """ 
    Return a copy of the configuration where an electron has been added to the
    given sub-shell (valence sub-shell by default).
    """
    out=copy.deepcopy(self)
    tss=ss
    if tss is None: 
      tss=max(out.subshells(indices=True))
    else:
      if not isinstance(tss,int): tss=rlf_ssindex(tss)
    if tss not in out.subshells(indices=True):
      out.v_conf[tss]=1
    else:
      if out.v_conf[tss] < rlf_ssoccupancy(tss):
        out.v_conf[tss]+=1
      else:
        if ss is None:
          out.v_conf[tss+1]=1
        else:
          msg='specified sub-shell is full: %s'%ss
          raise ErrorRAP(msg)
    return out

# ------------------------------------------------------------------------------

  def difference(self,conf):
    """ 
    Return dictionary of sub-shells and difference of occupation numbers.
    The sign of the difference is according to self - conf.
    """
    sss=self.subshells(indices=True)
    osss=conf.subshells(indices=True)
    out={}
    for ss in sss:
      w=self.occupancy(ss)
      if ss not in osss:
        out[rlf_subshell(ss)]=w
      else:
        dw=w-conf.occupancy(ss)
        if dw != 0: out[rlf_subshell(ss)]=dw
    for ss in osss:
      if ss in sss: continue
      out[rlf_subshell(ss)]=-conf.occupancy(ss)
    return out

# ------------------------------------------------------------------------------

  def to_jconfig(self):
    """ Return list of possible J configurations. """
    clist=[[]]
    for i in self:
      (n,l)=rlf_ssdecode(i)
      w=self[i]

      if l == 0:
        jss1=rlf_jssencode(n,l,+1)
        t1=jss1+str(w)
        tout=[]
        for tconf in clist:
          tw=tconf[:]+[t1]
          tout.append(tw)
        clist=tout[:]
        continue

      jss1=rlf_jssencode(n,l,-1)
      jss2=rlf_jssencode(n,l,+1)
      s1max=min(w,2*l)
      s1min=max(0,w-2*l-2)

      tmp=[]
      for is1 in range(s1min,s1max+1): 
        t1=jss1+str(is1)
        t2=jss2+str(w-is1)
        tmp.append([t1,t2])

      tout=[]
      for tconf in clist:
        for tt in tmp: 
          tw=tconf[:]+tt
          tout.append(tw)
      clist=tout[:]

    out=[]
    for cc in clist:
      cstr=' '.join(cc)
      out.append(rap_jconfig(cstr))
    return out

# ------------------------------------------------------------------------------

  def allowed_terms(self,encoded=True):
    """ 
    Return list of allowed LS terms for configuration.  Return value is a
    dictionary where the keys are the term strings and the values hold the
    number of terms having the same designation.  Set encoded=False to have
    each keys be a tuple containing 2S+1, L, and parity.
    """
    termlist={}
    termlist[(1,0,0)]=1
    for i in self.v_conf:
      w=self.v_conf[i]
      (n,l)=rlf_ssdecode(i)
      tl=rlf_equiv_lsterms(l,w)
      termlist=rlf_lsterm_couple(termlist,tl)
    if not encoded: return termlist

    out={}
    for (g,l,p),n in termlist.iteritems():
      t=str(g)+SPECLETTERS[l].upper()
      out[t]=n
    return out

# ------------------------------------------------------------------------------

  def eissner(self):
    """ Return configuration label in Eissner notation. """
    out=''
    for ss in self.subshells(indices=True):
      x=rlf_ss2eissner(ss)
      q=50+self[ss]
      out+=str(q)+x
    q=int(out[0])-5
    if q == 0:
      out=out[1:]
    else:
      out[0]=str(q)
    return out

################################################################################
class rap_states_ca(rap_states):
  """ Store set of configuration-averaged states. """
  def __init__(self,filename=None,**args):
    self.clear()
    if filename: self.loadfile(filename,**args)

# ------------------------------------------------------------------------------

  def clear(self):
    """ Erase all data. """
    rap_states.clear(self)
    self.v_cup='CA'
    self.v_conf={}

# ------------------------------------------------------------------------------

  def clear_levels(self):
    """ Just erase level data. """
    rap_states.clear_levels(self)
    self.v_cup='CA'
    self.v_conf={}

# ------------------------------------------------------------------------------

  def loadfile(self,filename,**args):
    """ Load in level information from level data file. """
    ft=rlf_filetype(filename)
    if ft == 'x006':
      self.load_xstar(filename,**args)
    else:
      msg='do not know how to read file: %s'%filename
      raise ErrorRAP(msg)

# ------------------------------------------------------------------------------

  def reindex(self,remap=None):
    """ 
    Reassign level indices.  Indices are remapped according to dictionary,
    remap, old (key) -> new (value).  If remap is None, then reindex according
    to energy order.
    """
    imap=rap_states.reindex(self,remap)
    self.v_conf=dict([(imap[i],val) for (i,val) in self.v_conf.items()])

# ------------------------------------------------------------------------------

  def configurations(self,levlist=None,imap=False):
    """ Return list of all configurations. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self._mapme,levlist)
    out=[]
    for i in levlist:
      if self.v_conf[i] in out: continue
      out.append(self.v_conf[i])
    return out

# ------------------------------------------------------------------------------

  def configuration(self,idx,imap=False):
    """ Return configuration of level with given index. """
    tidx=self._mapme(idx,imap)
    return self.v_conf[tidx]

# ------------------------------------------------------------------------------

  def level(self,idx,imap=False):
    """ Return all information for single level (conf,par,E) """
    tidx=self._mapme(idx,imap)
    out=(self.v_conf[tidx],self.v_parity[tidx], self.v_energy[tidx])
    return out

# ------------------------------------------------------------------------------

  def label(self,idx,imap=False,spacer=None):
    """ Return level label; optional spacer can be used instead of space. """
    tidx=self._mapme(idx,imap)
    out=str(self.v_conf[tidx])
    if spacer is not None: out=out.replace(' ',spacer)
    return out

# ------------------------------------------------------------------------------

  def weight(self,idx,imap=False):
    """ Return weight of level. """
    tidx=self._mapme(idx,imap)
    return self.v_conf[tidx].weight()

# ------------------------------------------------------------------------------

  def add_level(self,conf,energy,idx=False):
    """ Add new level to list """
    if idx is not False and idx in self:
      msg='index already exists: %s'%idx
      raise ErrorRAP(msg)
    tidx=idx
    if tidx is False: 
      tidx=1
      if len(self) > 0: tidx=max(self.indices())+1
    tconf=conf
    if not isinstance(tconf,rap_config): tconf=rap_config(tconf)
    if self.num_electrons() is None:
      self.v_nelec=tconf.num_electrons()
    else:
      if tconf.num_electrons() != self.num_electrons():
        msg='inconsistent number of electrons (expected: %s); conf: %s'%\
            (self.num_electrons(),tconf)
        raise ErrorRAP(msg)
    if tconf in self.v_conf.values():
      msg='duplicate state found: %s'%tconf
      raise ErrorRAP(msg)
    self.v_conf[tidx]=tconf
    self.v_parity[tidx]=tconf.parity()
    self.v_energy[tidx]=energy

# ------------------------------------------------------------------------------

  def add_levels(self,toadd,idxstart=None):
    """ 
    Add set of levels to structure. 

    The input list, toadd, contains tuples with: energy (Ry), index, and
    configuration.  If the indexing is to be automated, then set index in the
    tuples to -1.
    """
    if len(toadd) == 0: return
    toadd.sort()
    idxcur=1
    if len(self) > 0: idxcur=max(self.indices())+1
    if idxstart is not None:
      if idxstart < idxcur: 
        msg='invalid starting index: %s'%idxstart
        raise ErrorRAP(msg)
      idxcur=idxstart
    tmp={}
    nechk=None
    for (en,idx,conf) in toadd:
      tidx=idx
      if tidx in self or tidx in tmp: 
        msg='duplicate index: %s'%tidx
        raise ErrorRAP(msg)
      if tidx < 0:
        tidx=idxcur
        idxcur+=1
      tconf=conf
      if not isinstance(conf,rap_config): tconf=rap_config(conf)
      if nechk is None: nechk=tconf.num_electrons()
      if tconf.num_electrons() != nechk:
        msg='inconsistent number of electrons (expected: %s); conf: %s'%\
            (nechk,tconf)
        raise ErrorRAP(msg)
      tmp[tidx]=(en,tconf)
    self.v_nelec=nechk
    for idx,(en,conf) in tmp.iteritems():
      if conf in self.v_conf.values():
        msg='duplicate state found: %s'%conf
        raise ErrorRAP(msg)
      self.v_conf[idx]=conf
      self.v_parity[idx]=conf.parity()
      self.v_energy[idx]=en

# ------------------------------------------------------------------------------

  def set_lmap(self,other):
    """ Map to other rap_states_ca instance. """
    if not isinstance(other,rap_states_ca):
      msg='argument must be of type rap_states_ca'
      raise ErrorRAP(msg)
    if self.num_electrons() != other.num_electrons():
      msg='inconsistent number of electrons; expected/have: %s, %s'%\
          (self.num_electrons(),other.num_electrons())
      raise ErrorRAP(msg)
    tmap={}
    for i in other:
      k=None
      for j in self:
        if self[j] == other[i]:
          k=j
          break
      tmap[i]=k
    self.v_lmap=tmap
    for j in other:
      if j not in self.v_lmap: self.v_lmap[j]=None
    self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items() if v is not None])
    for i in self:
      if i not in self.v_lmap_r: self.v_lmap_r[i]=None

# ------------------------------------------------------------------------------

  def cfilter(self,clist,subset=None,imap=False):
    """ Filter list of states based on configurations in clist. """
    if not isinstance(clist,list) and not isinstance(clist,tuple): clist=[clist]
    tclist=[]
    for c in clist:
      tc=c
      if not isinstance(tc,rap_config): tc=rap_config(tc)
      tclist.append(tc)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.configuration(i) in tclist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def load_xstar(self,filename):
    """ 
    Load level information from xstar type 006 file.
    """
    # column id's
    id_energy=7
    id_j2=8
    id_ip=10
    id_g=12
    id_l=13
    id_z=14
    id_ion=16
    id_id=15
    id_conf=17

    # read level list from file
    ion=None
    ne=None
    ldata=[]
    have=[]
    xid0=None
    fin=open(filename,'r')
    for line in fin:
      if line.strip() == '': continue
      data=line.split()

      confstr=data[id_conf-1]
      xid=int(data[id_id-1])
      ion=int(data[id_ion-1])
      g=abs(int(data[id_g-1]))
      l=int(data[id_l-1])
      j2=int(float(data[id_j2-1]))-1
      w=g*(2*l+1)
      if g == 1 or l == 0:
        cup='XX'
      elif w < j2+1:
        cup='CA'
        if 'CA' not in have: have.append('CA')
      elif w == j2+1:
        cup='LS'
        if 'LS' not in have: have.append('LS')
      else:
        j2min=abs(g-1-2*l)
        j2max=g-1+2*l
        okay=False
        for j2t in range(j2min,j2max+1,2):
          if j2t == j2: okay=True
        if okay:
          cup='IC'
          if 'IC' not in have: have.append('IC')
        else:
          cup='??'
      e=float(data[id_energy-1])
      if e == 0.0:
        ne=rap_config(confstr).num_electrons()
        zeff=int(data[id_z-1])
        ion_thresh=float(data[id_ip-1])
      chk=confstr.lower()
      if 'su' in chk or 'sp' in chk: continue 
      if 'co' in chk or 'cn' in chk: continue
      if 'fake' in chk: continue
      ldata.append([e/RYDBERG,cup,confstr,g,l,j2,xid])
    fin.close()
    if ne is None: raise RuntimeError('cannot find ground state')

    # id's for ldata
    id_e=0
    id_cup=1
    id_conf=2
    id_g=3
    id_l=4
    id_j2=5
    id_xid=6

    # get number of electrons from ground state and standardize configurations
    for i in range(len(ldata)):
      ldata[i][id_conf]=rap_config(ldata[i][id_conf],ne=ne)

    # get total weight of each configuration
    weight={}
    wtot={}
    cups={}
    for i in range(len(ldata)):
      conf=str(ldata[i][id_conf])
      if conf not in wtot: wtot[conf]=ldata[i][id_conf].weight()
      if conf not in cups: cups[conf]=[]
      if ldata[i][id_cup] not in cups[conf]: cups[conf].append(ldata[i][id_cup])
      if ldata[i][id_cup] in ['XX','LS']:
        if conf not in weight:
          weight[conf]=ldata[i][id_g]*(2*ldata[i][id_l]+1)
        else:
          weight[conf]+=ldata[i][id_g]*(2*ldata[i][id_l]+1)
      else:
        if conf not in weight:
          weight[conf]=(ldata[i][id_j2]+1)
        else:
          weight[conf]+=(ldata[i][id_j2]+1)

    # check if configurations have a common coupling scheme
    for c in cups.keys():
      tmp=cups[c][:]
      if 'XX' in tmp: tmp.remove('XX')
      if len(tmp) > 1:
        raise RuntimeError('mixed coupling scheme not supported')
      if tmp == [] and len(have) == 1: tmp=have[:]
      cups[c]=tmp[:]
      if cups[c] == []: cups[c]=['CA']

    # check that 'CA' configurations have all possible levels
    for c in weight.keys():
      if 'CA' not in cups[c]: continue
      elif weight[c] > wtot[c]:
        raise RuntimeError('configuration, '+str(c)+', has too many levels')

    # add levels having CA coupling
    toadd=[]
    for lev in ldata:
      idx=lev[id_xid]
      conf=lev[id_conf]
      sc=str(conf)
      if 'CA' not in cups[sc] and cups[sc] != ['XX']: continue
      energy=lev[id_e]
      toadd.append((energy,idx,conf))
    if len(toadd) == 0:
      msg='no CA data found in file, %s'%filename
      raise ErrorRAP(msg)
    self.set_nuclear_charge(zeff)
    self.set_ionization_threshold(ion_thresh/RYDBERG)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def write_xstar(self,filename,subset=None,precision=5):
    """ 
    Write xstar type 6 file. Must have ionization potential and nuclear charge
    specified.
    """
    if self.ionization_threshold() == None:
      raise ErrorRAP('must define ionization threshold')
    if self.nuclear_charge() == None:
      raise ErrorRAP('must define nuclear charge')
    indices=self.indices()
    if subset is not None: indices=subset

    zeff=self.nuclear_charge()
    ne=self.num_electrons()
    ionno=xstar_ionencode(zeff,ne)

    aform='%'+str(8+precision)+'.'+str(precision)+'E'
    iform='%6i'
    prefix=iform%6+iform%13+iform%0+iform%4+iform%6
    fout=open(filename,'w')
    for i in indices:
      conf=self.configuration(i)
      e=self.energy(i)*RYDBERG
      wgt=float(self.weight(i))
      ip=self.ionization_threshold()*RYDBERG
      np=-1.0
      if ip > 0.0: 
        z1=float(zeff-ne+1)
        de=abs(self.ionization_threshold()-self.energy(i))
        np=z1/sqrt(de)
      nval=0
      if conf is not None: nval=rlf_ssdecode(conf.valence_subshell())[0]
      g,l=0,-1
      label=self.label(i,spacer='.')
      llen=len(label)
      line=prefix+iform%llen+' '+aform%e+aform%wgt+aform%np+aform%ip+' '+\
                  iform%nval+iform%g+iform%l+iform%zeff+iform%i+\
                  iform%ionno+' '+label+' %'
      fout.write(line+'\n')
    fout.close()

################################################################################
class rap_mix:
  """ Stores the mixing matrix from an olg file """
  def __init__(self,filename=None,levs=None,**args):
    self.clear()
    if levs is not None: self.v_levs=levs          # pointer to level object
    if filename: self.load_olg(filename,**args)

# ------------------------------------------------------------------------------

  def clear(self):
    self.v_levs=None
    self.v_mixmat={}               # mixing matrix, index is symmetry tuple
    self.v_symmap={}               # store symmetry tuple for each state

# ------------------------------------------------------------------------------

  def set_levs(self,levs):
    """ Set rap_states object. """
    if not isinstance(levs,rap_states):
      msg='levs must be a rap_states instance (or child instance)'
      raise ErrorRAP(msg)
    self.v_levs=levs

# ------------------------------------------------------------------------------

  def levs(self):
    """ Return rap_states object. """
    return self.v_levs

# ------------------------------------------------------------------------------

  def num_levels(self):
    """ Return number of levels. """
    return len(self.v_levs)

# ------------------------------------------------------------------------------

  def loadfile(self,filename,**args):
    """ Load in mixing data from data file. """
    ft=rlf_filetype(filename)
    if ft == 'olg':
      self.load_olg(filename,**args)
    else:
      msg='do not know how to read file: %s'%filename
      raise ErrorRAP(msg)

# ------------------------------------------------------------------------------

  def symmetries(self,states=None,imap=False):
    """
    Return list of symmetries in olg file. 
    
    Each symmetry is denoted by a tuple containing either (2S+1,L,parity)
    for LS-coupling or (2J,parity) for LSJ-coupling.  If states is given,
    only a subset of symmetries are returned corresponding to the terms in
    the list.
    """
    if imap and states is not None:
      states=map(self.v_levs._mapme,states)
    if not states:
      syms=self.v_mixmat.keys()
    else:
      try:
        syms=[]
        for t in states:
          s=self.symmetry(t)
          if s not in syms: syms.append(s)
      except:
        return False
    syms.sort()
    return syms

# ------------------------------------------------------------------------------

  def symmetry(self,idx,imap=False):
    """ Return symmetry for state with index, idx. """
    tidx=self.v_levs._mapme(idx,imap)
    return self.v_symmap[tidx]

# ------------------------------------------------------------------------------

  def entry(self,sym,imap=False):
    """ Return mixing matrix dictionary entry for single symmetry. """
    if sym not in self.v_mixmat.keys(): return False
    out=self.v_mixmat[sym]
    if not imap: return out
    keys=sorted(out.keys())
    keys=map(self.v_levs.lmap_reverse,keys)
    skeys=sorted(keys)
    tout={}
    for key,val in out.iteritems():
      newkey=self.v_levs.lmap_reverse(key)
      tmp=[]
      for i in range(len(keys)):
        tmp.append(out[key][keys.index(skeys[i])])
      tout[newkey]=tmp
    return tout

# ------------------------------------------------------------------------------

  def coeff(self,idx,imap=False):
    """ Return mixing coefficient dictionary for single state. """
    tidx=self.v_levs._mapme(idx,imap)
    s=self.symmetry(tidx)
    if not s: return False
    keys=self.v_levs.symfilter(s)
    coeffs=self.v_mixmat[s][tidx][:]
    out={}
    for k in keys: 
      kidx=k
      if imap: kidx=self.v_levs.lmap_reverse(kidx)
      out[kidx]=coeffs.pop(0)
    return out

# ------------------------------------------------------------------------------

  def matrix(self,symlist=None,imap=True):
    """ 
    Return 2D list containing mixing matrix for symmetries in symlist.
    """
    # get list of terms having symmetries
    tsymlist=copy.deepcopy(symlist)
    if not symlist:
      tsymlist=self.symmetries()
      termlist=self.v_levs.indices()
    elif symlist in self.symmetries():
      tsymlist=[symlist]
      termlist=self.v_levs.symfilter(symlist)
    else:
      termlist=[]
      for sym in symlist: termlist.extend(self.v_levs.symfilter(sym))
    if imap is True:
      termlist=map(self.v_levs.lmap_reverse,termlist)
    termlist.sort()

    # prepare final matrix and dummy row
    trow=[0 for t in termlist]
    tmat=[trow for t in termlist]

    # construct matrix
    for s in tsymlist:
      sterms=sorted(self.v_levs.symfilter(s,imap=imap))
      srow=self.entry(s,imap=imap)
      for tr in sterms:
        ridx=termlist.index(tr)
        row=trow[:]
        for i in range(len(termlist)):
          tc=termlist[i]
          if tc in sterms:
            j=sterms.index(tc)
            row[i]=srow[tr][j]
          else:
            row[i]=0.0
        tmat[ridx]=row
    return tmat

# ------------------------------------------------------------------------------

  def load_olg(self,filename,cup=None):
    """ Read information from olg file """
    # first some useful functions
    def get_sym(line,cup):
      # get symmetry tuple from line (global variable)
      try:
        if cup == 'LS':
          sym=(int(line[4:7]),int(line[7:10]),int(line[10:13]))
        else:
          sym=(int(line[8:11]),int(line[17:20]))
      except:
        return False
      return sym

    if self.levs() == None: self.v_levs=rlf_loadstates(filename,cup)
    tcup=self.v_levs.coupling()
    if cup is not None and cup != tcup:
      msg='given coupling scheme inconsistent with levels object'
      raise ErrorRAP(msg)
    cup=tcup

    # try to open olg file
    if rlf_filetype(filename) != 'olg': return False
    try:
      fin=open(filename,'r')
    except:
      return False

    # find start of mixing matrices
    line=fin.readline()
    while 'MATRIX' not in line: 
      if line == '': return False
      line=fin.readline()
    istart=51
    if cup == 'LSJ':
      istart=55
      line=fin.readline()
      while 'MATRIX' not in line: 
        if line == '': return False
        line=fin.readline()
    line=fin.readline()
    line=fin.readline()

    # read in mixing matrices
    termsym=[]                   # store symmetry for each term
    oldsym=get_sym(line,cup)
    idx=1
    termsym.append((idx,oldsym))
    cstr=line[istart:]
    alldata={}
    tmat=[]
    tdict={}
    done=False
    while not done:
      line=fin.readline()
      if 'CORE' in line: done=True
      sym=get_sym(line,cup)
      if not sym and not done:
        cstr+=line[istart:]
        continue
      cstr=cstr.strip()[:-8]
      coeffs=cstr.split()
      coeffs=map(float,coeffs)
      tdict[idx]=coeffs
      if done or sym != oldsym:
        nc=999999
        for t in tdict.keys(): nc=min(nc,len(tdict[t]))
        for t in tdict.keys(): tdict[t]=tdict[t][:nc]
        alldata[oldsym]=tdict
        tdict={}
      if sym:
        idx+=1
        cstr=line[istart:]
        oldsym=sym
        termsym.append((idx,oldsym))

    # get mapping from energy-order to symmetry-order
    while 'K*CM' not in line:
      if line == '': return False
      line=fin.readline()
    map_es=[]
    while True:
      line=fin.readline()
      tmp=line.split()
      try:
        it=int(tmp[1])
      except:
        break
      map_es.append(it)
    fin.close()

    # get reverse map
    map_se=[]
    for i in range(1,len(map_es)+1): 
      map_se.append(map_es.index(i)+1)

    # get list of symmetries for each term
    ttermsym=[]
    for (i,s) in termsym: ttermsym.append((map_se[i-1],s))
    ttermsym.sort()
    termsym={}
    for (i,s) in ttermsym: termsym[i]=s
    self.v_symmap=termsym

    # convert symmetry-ordered indices in mixing matrix to energy-ordered
    for s in alldata.keys():
      new={}
      oldterms=alldata[s].keys()
      oldterms.sort()
      newterms=[]
      for ot in oldterms: newterms.append(map_se[ot-1])
      newterms_sort=newterms[:]
      newterms_sort.sort()
      tmap=[]
      for nt in newterms: tmap.append(newterms_sort.index(nt)+1)
      for t in oldterms:
        mt=map_se[t-1]
        row=alldata[s][t]
        newrow=[0 for i in range(len(row))]
        for i in range(1,len(row)+1):
          newrow[tmap[i-1]-1]=row[i-1]
        new[mt]=newrow
      alldata[s]=new

    self.v_mixmat=alldata
    return True

################################################################################
class rap_states_lsj(rap_states):
  """ Store set of LSJ states. """
  def __init__(self,filename=None,**args):
    self.clear()
    if filename: self.loadfile(filename,**args)

# ------------------------------------------------------------------------------

  def _reprocess(self):
    """ Redefine senority-like tags. """
    if len(self) == 0: return
    indices=zip(*sorted(zip(self.v_energy.values(),self.v_energy.keys())))[1]
    allowedterms={}

    have={}
    tmp={}
    for idx in indices:
      (conf,term,sen,j2,par,en)=self.level(idx)

      # check for legal term
      if conf not in allowedterms: allowedterms[conf]=conf.allowed_terms()
      if term not in allowedterms[conf]:
        msg='invalid term for given configuration: %s %s'%(conf,term)
        raise ErrorRAP(msg)

      # check for legal level
      (g,l)=rlf_termdecode(term)
      if j2 < abs(2*l-g+1) or j2 > 2*l+g-1 or j2%2 != (2*l+g-1)%2:
        msg='invalid 2J value for term: %s %s 2J=%s'%(conf,term,j2)
        raise ErrorRAP(msg)

      # get and store seniority
      key=(conf,term)
      if key not in have: have[key]=0
      new=True
      for ij in range(1,have[key]+2):
        keyj=(conf,term,ij)
        if keyj not in tmp: break
        if j2 not in tmp[keyj]: 
          new=False
          break
      if new: 
        have[key]+=1
        tmp[keyj]=[]
      tmp[keyj].append(j2)
      self.v_seniority[idx]=ij
      if have[key] > allowedterms[conf][term]:
        msg='too many terms (%s) in configuration: %s %s; should have %s'%\
            (have[key],conf,term,allowedterms[conf][term])
        raise ErrorRAP(msg)

# ------------------------------------------------------------------------------

  def _parse_nist_line(self,line):
    """ 
    Function to extract level information (if any) from a NIST data file 
    (from webpage).  
    """
    if 'LIMIT' in line.upper():
      tmp=re.findall(r'(\d+\.\d+)',line)
      if not tmp: return False
      sym=line.split()[0].capitalize()
      zeff=None
      if sym in ATOMSYMBOLS: zeff=ATOMSYMBOLS.index(sym)
      ionth=float(tmp[-1])
      return (zeff,ionth)
    parts=line.split('|')
    if parts[0] == line: return False
    if parts[0].strip() == '':
      conf=None                  # None means no conf info
    else:
      tmp=re.findall(r'(\d+[a-z]\d*)',parts[0])
      if tmp == []: 
        conf=False               # False means can't understand conf
      else:
        conf=''
        for t in tmp: 
          l=re.findall(r'([a-z])',t)[0]
          l=SPECLETTERS.index(l.upper())
          tmp=re.findall(r'(\d+)',t)
          if len(tmp) == 1: t=t+'1'
          conf=conf+t+' '
        conf=conf[:-1]
    if parts[1].strip() == '':
      term=None                  # None means no term info
    else:
      tmp=re.findall(r'(\d+[A-Z])',parts[1])
      if tmp == []:
        term=False               # False means can't understand term
      else:
        term=tmp[0]
    tmp=re.findall(r'(\*)',parts[1])
    tmp=re.findall(r'(\d+)/2',parts[2])
    if tmp == []:
      tmp=re.findall(r'(\d+)',parts[2])
      if tmp == []: return False
      jay=float(tmp[0])
    else:
      jay=0.5*float(tmp[0])
    tmp=re.findall(r'(\d+\.\d+)',parts[3])
    energy=None
    if tmp != []: energy=float(tmp[0].replace('+x',''))
    j2=int(2.0*jay)
    return (conf,term,j2,energy)

# ------------------------------------------------------------------------------

  def clear(self):
    """ Erase all data. """
    rap_states.clear(self)
    self.v_cup='LSJ'
    self.v_conf={}
    self.v_term={}
    self.v_seniority={}
    self.v_j2={}

# ------------------------------------------------------------------------------

  def clear_levels(self):
    """ Just erase level data. """
    rap_states.clear_levels(self)
    self.v_cup='LSJ'
    self.v_conf={}
    self.v_term={}
    self.v_seniority={}
    self.v_j2={}

# ------------------------------------------------------------------------------

  def loadfile(self,filename,**args):
    """ Load in level information from level data file. """
    ft=rlf_filetype(filename)
    if ft == 'adf04':
      self.load_adf04(filename,**args)
    elif ft == 'olg':
      self.load_olg(filename,**args)
    elif ft == 'nist':
      self.load_nist(filename,**args)
    elif ft == 'chianti_elvlc':
      self.load_chianti(filename,**args)
    elif ft == 'x006':
      self.load_xstar(filename,**args)
    elif ft == 'atomdb_ev':
      self.load_atomdb(filename,**args)
    else:
      msg='do not know how to read file: %s'%filename
      raise ErrorRAP(msg)

# ------------------------------------------------------------------------------

  def reindex(self,remap=None):
    """ 
    Reassign level indices.  Indices are remapped according to dictionary,
    remap, old (key) -> new (value).  If remap is None, then reindex according
    to energy order.
    """
    imap=rap_states.reindex(self,remap)
    self.v_conf=dict([(imap[i],val) for (i,val) in self.v_conf.items()])
    self.v_term=dict([(imap[i],val) for (i,val) in self.v_term.items()])
    self.v_seniority=dict([(imap[i],val) for (i,val) in \
                          self.v_seniority.items()])
    self.v_j2=dict([(imap[i],val) for (i,val) in self.v_j2.items()])

# ------------------------------------------------------------------------------

  def configurations(self,levlist=None,imap=False):
    """ Return list of all configurations. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self._mapme,levlist)
    out=[]
    for i in levlist:
      if self.v_conf[i] in out: continue
      out.append(self.v_conf[i])
    return out

# ------------------------------------------------------------------------------

  def configuration(self,idx,imap=False):
    """ Return configuration of level with given index. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_conf[tidx]

# ------------------------------------------------------------------------------

  def terms(self,levlist=None):
    """ Return list of all terms. """
    if levlist is None: levlist=self.indices()
    out=[]
    for i in levlist:
      if self.v_term[i] in out: continue
      out.append(self.v_term[i])
    return out

# ------------------------------------------------------------------------------

  def term(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_term[tidx]

# ------------------------------------------------------------------------------

  def seniority(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_seniority[tidx]

# ------------------------------------------------------------------------------

  def j2range(self,levlist=None,imap=False):
    """ Return tuple containing the minimum and maximum 2J values. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self.mapme,levlist)
    j2s=[self.v_j2[i] for i in levlist if i is not None]
    return (min(j2s),max(j2s))

# ------------------------------------------------------------------------------

  def j2(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_j2[tidx]

# ------------------------------------------------------------------------------

  def jrange(self,levlist=None,imap=False):
    """ Return tuple containing the minimum and maximum 2J values. """
    (j2min,j2max)=self.j2range(levlist,imap)
    return (0.5*float(j2min),0.5*float(j2max))

# ------------------------------------------------------------------------------

  def j(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return 0.5*float(self.v_j2[tidx])

# ------------------------------------------------------------------------------

  def level(self,idx,imap=False):
    """ Return all information for single level (conf,term,alpha,j2,par,E) """
    tidx=self._mapme(idx,imap)
    if tidx is None: return (None,None,None,None,None,None)
    out=(self.v_conf[tidx],self.v_term[tidx],self.v_seniority[tidx],\
         self.v_j2[tidx],self.v_parity[tidx], self.v_energy[tidx])
    return out

# ------------------------------------------------------------------------------

  def label(self,idx,imap=False,spacer=None):
    """ Return level label; optional spacer can be used instead of space. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    j2=self.j2(idx,imap=imap)
    if j2%2 == 0:
      jstr='_'+str(j2/2)
    else:
      jstr='_'+str(j2)+'/2'
    out=str(self.v_conf[tidx])+' '+self.v_term[tidx]+jstr
    if self.v_seniority[tidx] > 1: out+='#'+str(self.v_seniority[tidx])
    if spacer is not None: out=out.replace(' ',spacer)
    return out

# ------------------------------------------------------------------------------

  def weight(self,idx,imap=False):
    """ Return weight of level. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.j2(idx,imap=imap)+1

# ------------------------------------------------------------------------------

  def add_level(self,conf,term,j2,energy,idx=False):
    """ Add new level to list """
    if idx is not False and idx in self:
      msg='index already exists: %s'%idx
      raise ErrorRAP(msg)
    tidx=idx
    if tidx is False: 
      tidx=1
      if len(self) > 0: tidx=max(self.indices())+1
    tconf=conf
    if not isinstance(tconf,rap_config): tconf=rap_config(tconf)
    if self.num_electrons() is None:
      self.v_nelec=tconf.num_electrons()
    else:
      if tconf.num_electrons() != self.num_electrons():
        msg='inconsistent number of electrons; expect/have: %s, %s'%\
            (self.num_electrons(),tconf.num_electrons())
        raise ErrorRAP(msg)
    (g,l)=rlf_termdecode(term)     # will raise error if bad term string
    if j2 < abs(2*l-g+1) or j2 > 2*l+g-1 or j2%2 != (2*l+g-1)%2:
      msg='invalid 2J value for term: %s 2J=%s'%(term,j2)
      raise ErrorRAP(msg)
    self.v_conf[tidx]=tconf
    self.v_term[tidx]=term
    self.v_seniority[tidx]=-1
    self.v_j2[tidx]=j2
    self.v_parity[tidx]=tconf.parity()
    self.v_energy[tidx]=energy
    self._reprocess()

# ------------------------------------------------------------------------------

  def add_levels(self,toadd,idxstart=None):
    """ 
    Add set of levels to structure. 

    The input list, toadd, contains tuples with: energy (Ry), index,
    configuration, term, and 2J.  If the indexing is to be automated, then set 
    index in the tuples to -1.
    """
    if len(toadd) == 0: return
    toadd.sort()
    idxcur=1
    if len(self) > 0: idxcur=max(self.indices())+1
    if idxstart is not None:
      if idxstart < idxcur: 
        msg='invalid starting index: %s'%idxstart
        raise ErrorRAP(msg)
      idxcur=idxstart
    tmp={}
    nechk=None
    for (en,idx,conf,term,j2) in toadd:
      tidx=idx
      if tidx in self or tidx in tmp: 
        msg='duplicate index: %s'%tidx
        raise ErrorRAP(msg)
      if tidx < 0:
        tidx=idxcur
        idxcur+=1
      tconf=conf
      if not isinstance(conf,rap_config): tconf=rap_config(conf)
      if nechk is None: nechk=tconf.num_electrons()
      if tconf.num_electrons() != nechk:
        msg='inconsistent number of electrons; expected/have: %s, %s'%\
            (nechk,tconf.num_electrons())
        raise ErrorRAP(msg)
      try:                       # raise error if bad term string
        (g,l)=rlf_termdecode(term)
      except ErrorRAP,msg:
        msg='problem with %s %s 2J=%s... %s'%(conf,term,j2,msg)
        warnings.warn(msg)
        continue
      if j2 not in range(abs(2*l-g+1),2*l+g,2):
        msg='invalid 2J value for term: %s %s 2J=%s'%(conf,term,j2)
        raise ErrorRAP(msg)
      tmp[tidx]=(en,tconf,term,j2)
    self.v_nelec=nechk
    for idx,(en,conf,term,j2) in tmp.iteritems():
      self.v_conf[idx]=conf
      self.v_term[idx]=term
      self.v_seniority[idx]=-1
      self.v_j2[idx]=j2
      self.v_parity[idx]=conf.parity()
      self.v_energy[idx]=en
    self._reprocess()

# ------------------------------------------------------------------------------

  def set_lmap(self,other,simple=False,minimize=False):
    """ Map to other rap_states_lsj instance. """
    if not isinstance(other,rap_states_lsj):
      msg='argument must be of type rap_states_lsj'
      raise ErrorRAP(msg)
    if self.num_electrons() != other.num_electrons():
      msg='inconsistent number of electrons; expected/have: %s, %s'%\
          (self.num_electrons(),other.num_electrons())
      raise ErrorRAP(msg)

    if simple:
      tmap={}
      for i in other:
        k=None
        for j in self:
          if self[j] == other[i]:
            k=j
            break
        tmap[i]=k
      self.v_lmap=tmap
      self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items()])
      return

    # group level indices by (conf,term,2J); seniority numbers not trusted
    fromlist={}
    for i in self:
      key=(self.configuration(i),self.term(i),self.j2(i))
      if key not in fromlist: fromlist[key]=[]
      fromlist[key].append(i)
    tolist={}
    for j in other:
      key=(other.configuration(j),other.term(j),other.j2(j))
      if key not in tolist: tolist[key]=[]
      tolist[key].append(j)

    if minimize:
      self.v_lmap=self._map_by_energy(other,fromlist,tolist)
    else:
      self.v_lmap=self._map_straight(other,fromlist,tolist)
    for j in other:
      if j not in self.v_lmap: self.v_lmap[j]=None
    self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items() if v is not None])
    for i in self:
      if i not in self.v_lmap_r: self.v_lmap_r[i]=None

# ------------------------------------------------------------------------------

  def cfilter(self,clist,subset=None,imap=False):
    """ Filter list of states based on configurations in clist. """
    if not isinstance(clist,list) and not isinstance(clist,tuple): clist=[clist]
    tclist=[]
    for c in clist:
      tc=c
      if not isinstance(tc,rap_config): tc=rap_config(tc)
      tclist.append(tc)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.configuration(i) in tclist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def tfilter(self,tlist,alpha=None,subset=None,imap=False):
    """ 
    Filter list of states based on terms in tlist; allow for optional 
    input for the seniority (alpha).
    """
    if not isinstance(tlist,list) and not isinstance(tlist,tuple): tlist=[tlist]
    if alpha is not None:
      if not isinstance(alpha,list) and not isinstance(alpha,tuple): 
        alpha=[alpha]

    ttlist=[]
    for t in tlist:
      tt=t
      if isinstance(tt,tuple): tt=rlf_termencode(tt[0],tt[1])
      ttlist.append(tt)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if alpha is not None and self.seniority(i) not in alpha: continue
      if self.term(i) in ttlist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def j2filter(self,j2list,subset=None,imap=False):
    """ Filter list of states based on 2J values in j2list. """
    if not isinstance(j2list,list) and not isinstance(j2list,tuple): 
      j2list=[j2list]
    tj2list=map(int,j2list)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.j2(i) in tj2list: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def jfilter(self,jlist,subset=None,imap=False):
    """ Filter list of states based on J values in j2list. """
    if not isinstance(jlist,list) and not isinstance(jlist,tuple): jlist=[jlist]
    tjlist=map(float,jlist)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.j(i) in tjlist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def symfilter(self,j2,p=None,subset=None,imap=False):
    """ Filter list of states based on symmetry. """
    if p is None:
      try:
        (jj2,pp)=j2
      except:
        msg='do not understand symmetry'
        raise ErrorRAP(msg)
    else:
      jj2=j2
      pp=p

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      tj2=self.j2(i)
      tp=self.parity(i)
      if tj2 != jj2 or tp != pp: continue
      out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def LS_average(self,omit_incomplete=True):
    """ 
    Average levels into LS terms contained in rap_states_ls object. 

    Indices will be set to energy-ordered.
    """
    agg={}
    for i in self:
      key=(self.v_conf[i],self.v_term[i],self.v_seniority[i])
      if key not in agg: agg[key]=(0.0,0.0)
      (esum,wgt)=agg[key]
      esum+=float(self.v_j2[i]+1)*self.v_energy[i]
      wgt+=float(self.v_j2[i]+1)
      agg[key]=(esum,wgt)

    emin=min([esum/wgt for (esum,wgt) in agg.values()])
    toadd=[]
    for (conf,term,alfa),(esum,wgt) in agg.iteritems():
      if omit_incomplete and rlf_termweight(term) > int(wgt): continue
      e=esum/wgt
      toadd.append((e-emin,-1,conf,term))

    out=rap_states_ls()
    out.set_nuclear_charge(self.nuclear_charge())
    out.set_core_energy(self.core_energy())
    out.set_ionization_threshold(self.ionization_threshold())
    out.add_levels(toadd)
    return out

# ------------------------------------------------------------------------------

  def search(self,conf=None,term=None,j2=None,parity=None,alpha=None,\
             energy=None):
    """ 
    Look for matching levels in structure based on configuration, term, and
    2J values.  If duplicate terms are present and alpha is not given, a list
    of all matching indices are returned.  If alpha is given, the index matching
    that seniority is returned.  If energy is given (and alpha is not), then
    the senority closest to that energy is returned.  If conf and term are 
    both None, then return matching index based on j2, parity, and energy order
    given by alpha.
    """
    if j2 is None: 
      if conf is None or term is None:
        msg="if no 2J is specified, must provide configuration and term"
        raise ErrorRAP(msg)
      if alpha is None: alpha=1
      levs=self.cfilter(conf)
      levs=self.tfilter(term,alpha,subset=levs)
      return levs

    levs=self.j2filter(j2)
    if conf is None or term is None:
      if conf is not None: parity=conf.parity()
      if parity is None: 
        raise RuntimeError('need parity when conf or term is None')
      levs=self.pfilter(parity,subset=levs)
      if levs == []: raise RuntimeError('matching level not found')
      if alpha is not None:
        return levs[alpha-1]
      elif energy is not None:
        tmp=[]
        for i in levs: tmp.append((abs(energy-self.energy(i)),i))
        return min(tmp)[1]
      return levs

    levs=self.cfilter(conf,subset=levs)
    levs=self.tfilter(term,subset=levs)
    if levs == []: raise RuntimeError('matching level not found')
    if len(levs) == 1: return levs[0]
    if alpha is not None:
      for i in levs: 
        if self.seniority(i) == alpha: return i
      raise RuntimeError('matching level not found')
    elif energy is not None:
      tmp=[]
      for i in levs: tmp.append((abs(energy-self.energy(i)),i))
      return min(tmp)[1]
    return levs

# ------------------------------------------------------------------------------

  def load_olg(self,filename):
    """ Load in level information from an AUTOSTRUCTURE olg file """
    fin=open(filename,'r')

    # first get closed subshells and number of configurations
    closed=[]
    while 1:
      line=fin.readline()
      if line == '': 
        msg='configuration list not found'
        raise ErrorRAP(msg)
      if 'CONFIGURATIONS' not in line: continue
      tmp=re.findall(r'IN(.+)CONFIGURATIONS',line)
      if tmp: break
    nconf=int(tmp[0])
    while 1:
      tmp=re.findall(r'(\d+)\s+(\d+)',line)
      if tmp == []: break
      closed.extend(tmp)
      line=fin.readline()

    # get list of all open subshells
    while 1:
      line=fin.readline()
      if line == '': 
        msg='could not find sub-shell list'
        raise ErrorRAP(msg)
      if 'INTERNAL ORBITALS' in line: break
    torbs=[]
    while 1:
      line=fin.readline()
      if 'CONFIGURATION' in line: break
      if 'N L' in line: 
        tmp=re.findall(r'(\d+)\s+(\d+)',line)
        if tmp: torbs.extend(tmp)
    orbs=[]
    for o in torbs:
      if o not in closed: orbs.append(o)
    base={}
    for o in orbs: base[o]=0

    ne=None
    conflist=[]
    for iconf in range(nconf):
      while 1:
        if line == '': 
          msg='configuration list not found'
          raise ErrorRAP(msg)
        if 'COMBINATIONS' in line:
          tmp=re.findall(r'COMBINATIONS(.+)',line)
          if tmp: break
        line=fin.readline()
      tmp=tmp[0]
      while 1:
        line=fin.readline()
        if 'SLATER' in line: break
        tmp+=line[:-1]
      tconf=base.copy()
      nl=re.findall(r'(\d+)\s+(\d+)',tmp)
      if not nl: 
        msg='could not understand %s configuration'%rlf_ordinal(iconf)
        raise ErrorRAP(msg)
      for ss in nl: tconf[ss]+=1
      aconf=''
      for o in orbs:
        ss=str(o[0])+SPECLETTERS[int(o[1])].lower()
        n=tconf[o]
        aconf+=ss+str(n)+' '
      conflist.append(rap_config(aconf,ne=ne))
      if ne is None: ne=conflist[0].num_electrons()

    # get nuclear charge
    while 1:
      if line == '':
        msg='file ended before retrieval of nuclear charge'
        raise ErrorRAP(msg)
      line=fin.readline()
      if 'ATOMIC NUMBER' in line: break
    tmp=re.findall(r'ATOMIC NUMBER\s*(\d+)',line)
    zeff=int(tmp[0])

    # get level information
    while 1:
      line=fin.readline()
      if line == '': 
        msg='no level data found'
        raise ErrorRAP(msg)
      if '(EK-E1)/RY' in line: break
    tmp=re.findall(r'(-?\d+\.\d+)\s*$',line)
    core_energy=None
    if tmp: core_energy=float(tmp[0])
    toadd=[]
    while 1:
      line=fin.readline()
      tmp=line.split()
      if tmp == []: break
      try:
        idx=int(tmp[0])
      except:
        break
      term=rlf_termencode(abs(int(tmp[4])),int(tmp[5]))
      j2=int(tmp[6])
      iconf=int(tmp[7])-1
      if iconf < 0: continue
      conf=conflist[iconf]
      energy=float(tmp[8])
      toadd.append((energy,idx,conf,term,j2))
    fin.close()
    self.set_nuclear_charge(zeff)
    self.set_core_energy(core_energy)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_adf04(self,filename):
    """ Load in level information from an ADAS adf04 file """
    fin=open(filename,'r')

    # basic info
    line=fin.readline()
    tmp=re.findall(r'([\d.]+)',line)
    z=int(tmp[1])
    q=int(tmp[2])
    ip=float(tmp[3])
    ne=z-q+1

    # read in level data
    isic=True
    isls=True
    isca=True
    toadd=[]
    while 1:
      line=fin.readline()
      if '-1' in line: break
      tmp=re.findall(r'(\d+)(.+)\((.+)\)(.+)\((.+)\)\s+(\d+[\.\d]*)',line)[0]
      idx=int(tmp[0])
      conf=rap_config(tmp[1],ne=ne)
      g=int(tmp[2])
      l=int(tmp[3])
      j2=int(2.*float(tmp[4]))
      e=float(tmp[5])/INVERSECM
      toadd.append((e,idx,conf,rlf_termencode(g,l),j2))

      cwgt=conf.weight()
      lwgt=g*(2*l+1)
      wgt=j2+1
      if wgt != cwgt: isca=False
      if wgt != lwgt: isls=False
      if wgt not in  range(abs(2*l-g+1)+1,2*l+g+1,2): isic=False
    fin.close()

    # check for correct coupling scheme
    if not isic:
      if isca: 
        msg='adf04 file seems to be configuration averaged'
        raise ErrorRAP(msg)
      elif isls:
        msg='adf04 file seems to be LS averaged'
        raise ErrorRAP(msg)
      else:
        msg='adf04 file not LSJ or a mistake in the J-values'
        raise ErrorRAP(msg)

    self.set_nuclear_charge(z)
    self.set_ionization_threshold(ip)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_nist(self,filename):
    """ Load in NIST data file (text file taken from website). """
    ne=None
    zeff=None
    ip=None
    toadd=[]
    ionth=[]
    fin=open(filename,'r')
    skipped=0
    for line in fin:
      tmp=self._parse_nist_line(line)
      if tmp:
        if len(tmp) == 2:
          if tmp[0] is not None: zeff=tmp[0]
          ionth.append(tmp[1])
        else:
          j2=tmp[2]
          en=tmp[3]
          if tmp[0] is False:
            conf=None
            term=None
            skipped+=1
            continue
          elif tmp[0] is None:
            if conf is None or tmp[1] is False: 
              conf=None
              term=None
              skipped+=1
              continue
            conf=toadd[-1][2]
            term=toadd[-1][3]
          else:
            if tmp[1] is False: 
              conf=None
              term=None
              skipped+=1
              continue
            conf=rap_config(tmp[0],ne=ne)
            term=tmp[1]
          if en == 0.0: ne=conf.num_electrons()
          toadd.append((en,-1,conf,term,j2))
    fin.close()

    # remove states with no energy
    for i in range(len(toadd)-1,-1,-1):
      if toadd[i][0] is None: toadd.pop(i)

    # any levels skipped
    if skipped > 0:
      plural=''
      if skipped > 1: plural='s'
      msg='%s level%s not understood and skipped'%(skipped,plural)
      warnings.warn(msg)

    if len(ionth) > 0: ip=min(ionth)
    self.set_nuclear_charge(zeff)
    self.set_ionization_threshold(ip)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_chianti(self,filename,priority="experimental",only=False):
    """ 
    Read energy levels from CHIANTI file.  CHIANTI energy level files contain
    both theoretical and experimental energies.  The priority option allows
    you to select which energy takes precedence (if available).  The only flag
    will restrict this function to storing only levels where the priority
    energy is present.
    """
    id_idx=0         # level index
    id_k=1           # configuration index
    id_conf=2        # configuration string
    id_g=3           # 2S+1
    id_lval=4        # L (integer)
    id_l=5           # L (string)
    id_j=6           # J
    id_wgt=7         # weight (2J+1)
    id_obcm=8        # observed energy (1/cm)
    id_obry=9        # observed energy (Ry); bad conversion, use 1/cm value
    id_thcm=10       # theoretical energy (1/cm)
    id_thry=11       # theoretical energy (Ry); bad conversion, use 1/cm value

    # attempt to get nuclear charge from filename
    tmp=re.findall(r'(\w+)_\d+\.elvlc',filename)
    if tmp:  
      zeff=ATOMSYMBOLS.index(tmp[0].capitalize())
      self.set_nuclear_charge(zeff)

    ne=None
    toadd=[]
    fin=open(filename,'r')
    for line in fin:
      tmp=line.split()
      if tmp[0] == '-1': break
      if len(tmp) != 12:
        raise RuntimeError('cannot undertand CHIANTI energy level file')
      idx=int(tmp[id_idx])
      conf=rap_config(tmp[id_conf],ne=ne)
      parity=conf.parity()
      g=int(tmp[id_g])
      l=int(tmp[id_lval])
      term=rlf_termencode(g,l)
      j2=int(tmp[id_wgt])-1
      ob=float(tmp[id_obcm])/INVERSECM
      th=float(tmp[id_thcm])/INVERSECM
      if priority[0] == 't':
        en=th
        if en == 0.0 and idx > 1: 
          if only is True: continue
          en=ob
      else:
        en=ob
        if en == 0.0 and idx > 1: 
          if only is True: continue
          en=th
      if ne is None: ne=conf.num_electrons()
      toadd.append((en,idx,conf,term,j2))
    fin.close()
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_xstar(self,filename):
    """ 
    Load level information from xstar type 006 file.
    """
    # column id's
    id_energy=7
    id_j2=8
    id_ip=10
    id_g=12
    id_l=13
    id_z=14
    id_ion=16
    id_id=15
    id_conf=17

    # read level list from file
    ion=None
    ne=None
    ldata=[]
    have=[]
    xid0=None
    fin=open(filename,'r')
    for line in fin:
      if line.strip() == '': continue
      data=line.split()

      confstr=data[id_conf-1]
      xid=int(data[id_id-1])
      ion=int(data[id_ion-1])
      g=abs(int(data[id_g-1]))
      l=int(data[id_l-1])
      j2=int(float(data[id_j2-1]))-1
      w=g*(2*l+1)
      if g == 1 or l == 0:
        cup='XX'
      elif w < j2+1:
        cup='CA'
        if 'CA' not in have: have.append('CA')
      elif w == j2+1:
        cup='LS'
        if 'LS' not in have: have.append('LS')
      else:
        j2min=abs(g-1-2*l)
        j2max=g-1+2*l
        okay=False
        for j2t in range(j2min,j2max+1,2):
          if j2t == j2: okay=True
        if okay:
          cup='IC'
          if 'IC' not in have: have.append('IC')
        else:
          cup='??'
      e=float(data[id_energy-1])
      if e == 0.0:
        ne=rap_config(confstr).num_electrons()
        zeff=int(data[id_z-1])
        ion_thresh=float(data[id_ip-1])
      chk=confstr.lower()
      if 'su' in chk or 'sp' in chk: continue 
      if 'co' in chk or 'cn' in chk: continue
      if 'fake' in chk: continue
      ldata.append([e/RYDBERG,cup,confstr,g,l,j2,xid])
    fin.close()
    if ne is None: raise RuntimeError('cannot find ground state')

    # id's for ldata
    id_e=0
    id_cup=1
    id_conf=2
    id_g=3
    id_l=4
    id_j2=5
    id_xid=6

    # get number of electrons from ground state and standardize configurations
    for i in range(len(ldata)):
      ldata[i][id_conf]=rap_config(ldata[i][id_conf],ne=ne)

    # get total weight of each configuration
    weight={}
    wtot={}
    cups={}
    for i in range(len(ldata)):
      conf=str(ldata[i][id_conf])
      if conf not in wtot: wtot[conf]=ldata[i][id_conf].weight()
      if conf not in cups: cups[conf]=[]
      if ldata[i][id_cup] not in cups[conf]: cups[conf].append(ldata[i][id_cup])
      if ldata[i][id_cup] in ['XX','LS']:
        if conf not in weight:
          weight[conf]=ldata[i][id_g]*(2*ldata[i][id_l]+1)
        else:
          weight[conf]+=ldata[i][id_g]*(2*ldata[i][id_l]+1)
      else:
        if conf not in weight:
          weight[conf]=(ldata[i][id_j2]+1)
        else:
          weight[conf]+=(ldata[i][id_j2]+1)

    # check if configurations have a common coupling scheme
    for c in cups.keys():
      tmp=cups[c][:]
      if 'XX' in tmp: tmp.remove('XX')
      if len(tmp) > 1:
        raise RuntimeError('mixed coupling scheme not supported')
      if tmp == [] and len(have) == 1: tmp=have[:]
      cups[c]=tmp[:]
      if cups[c] == []: cups[c]=['IC']

    # check that 'IC' configurations have all possible levels
    for c in weight.keys():
      if 'IC' not in cups[c]: continue
      elif weight[c] > wtot[c]:
        raise RuntimeError('configuration, '+str(c)+', has too many levels')

    # add levels having IC coupling
    toadd=[]
    for lev in ldata:
      idx=lev[id_xid]
      conf=lev[id_conf]
      sc=str(conf)
      if 'IC' not in cups[sc] and cups[sc] != ['XX']: continue
      term=rlf_termencode(lev[id_g],lev[id_l])
      j2=lev[id_j2]
      energy=lev[id_e]
      toadd.append((energy,idx,conf,term,j2))
    if len(toadd) == 0:
      msg='no LSJ data found in file, %s'%filename
      raise ErrorRAP(msg)
    self.set_nuclear_charge(zeff)
    self.set_ionization_threshold(ion_thresh/RYDBERG)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_atomdb(self,filename):
    """ 
    Load level information from ATOMDB FITS file.
    """
    # need pyfits
    rlf_have_pyfits(stop=True)

    a=pyfits.open(filename)
    adat=a[1].data
    header=a[1].header

    # check file type
    if header['HDUCLAS1'] != 'E_LEVEL':
      msg="invalid HDUCLAS1 in FITS file"
      raise ErrorRAP(msg)

    ionpot=None
    if 'ION_POT' in a[1].header:
      ionpot=float(a[1].header['ION_POT'])

    ionstr=header['ION_NAME']
    (z,ne)=rlf_iondecode(ionstr)
    q=z-ne

    # read in level data
    toadd=[]
    idx=0
    wgts={}                              # total configuration weight
    skipped=0
    for i in range(len(adat)):
      idx+=1
      g=int(2.*float(adat.field('S_QUAN')[i]))+1
      l=int(adat.field('L_QUAN')[i])
      j2=int(adat.field('LEV_DEG')[i])-1
      if g < 0:                                   # not LSJ
        skipped+=1
        continue
      e=float(adat.field('ENERGY')[i])            # in eV
      conf=rap_config(adat.field('ELEC_CONFIG')[i],ne=ne)

      if conf not in wgts: wgts[conf]=0
      wgts[conf]+=j2+1

      term=rlf_termencode(g,l)
      toadd.append((e/RYDBERG,idx,conf,term,j2))
    a.close()

    # check configuration weights
    for c in wgts:
      if c.weight() != wgts[c]:
        msg='incomplete configuration: %s'%str(c)
        warnings.warn(msg)

    if skipped > 0:
      msg='%d levels skipped; wrong coupling scheme'%skipped
      warnings.warn(msg)

    self.set_nuclear_charge(z)
    if ionpot is not None: self.set_ionization_threshold(ionpot)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def write_atomdb(self,filename,eref=None,comments=None):
    """ 
    Write level information to ATOMDB FITS file (from Adam Foster).
    """
    # need pyfits
    rlf_have_pyfits(stop=True)

    # prep comments
    if comments is None:
      comments=[]
    elif isinstance(comments,str):
      comments=[comments]
    else:
      comments=list(comments)

    # reference for energy levels
    if eref is None: eref=''

    # assemble energy level data
    elec_config=[]
    energy=[]
    e_error=[]
    n_quan=[]
    l_quan=[]
    s_quan=[]
    lev_deg=[]
    phot_type=[]
    phot_par=[]
    energy_ref=[]
    phot_ref=[]
    for ilev in self:
      (conf,term,alfa,j2,par,en)=self.level(ilev)
      valn=rlf_ssdecode(conf.valence_subshell())[0]
      (g,l)=rlf_termdecode(term)
      s=0.5*float(g-1)

      elec_config.append(str(conf))
      energy.append(RYDBERG*en)
      e_error.append(0.0)
      n_quan.append(valn)
      l_quan.append(l)
      s_quan.append(s)
      lev_deg.append(j2+1)
      phot_type.append(0)
      phot_par.append([0.0 for i in range(20)])
      energy_ref.append(eref)
      phot_ref.append('')

    # ion string
    z=self.nuclear_charge()
    ne=self.num_electrons()
    q=z-ne
    ionstr=rlf_ionencode(z,ne)

    # ionization potential
    ionpot=self.ionization_threshold()
    if ionpot is not None: ionpot*=RYDBERG

    #primary HDU, hdu0
    hdu0 = pyfits.PrimaryHDU()
    now = datetime.datetime.utcnow()

    hdu0.header.update('DATE', now.strftime('%d/%m/%y'))
    hdu0.header.update('E_LEVEL', "Atomic Energy Levels")
    hdu0.header.update('FILENAME', "Python routine")
    hdu0.header.update('ORIGIN', "ATOMDB",comment=os.environ['USER']+\
                       ", AtomDB project")
    hdu0.header.update('HDUCLASS', "ATOMIC",comment="Atomic Data")
    hdu0.header.update('HDUCLAS1', "E_LEVEL",comment="Atomic Energy levels")
    hdu0.header.update('HDUVERS', "1.1.0",comment="Version of datafile")

    #secondary HDU, hdu1:
    hdu1 = pyfits.new_table(pyfits.ColDefs(
    	      [pyfits.Column(name='ELEC_CONFIG',
    			     format='40A',
    			     array=elec_config),
    	       pyfits.Column(name='ENERGY',
    			     format='1E',
    			     unit='eV',
    			     array=energy),
    	       pyfits.Column(name='E_ERROR',
    			     format='1E',
    			     unit='eV',
    			     array=e_error),
    	       pyfits.Column(name='N_QUAN',
    			     format='1J',
    			     array=n_quan),
    	       pyfits.Column(name='L_QUAN',
    			     format='1J',
    			     array=l_quan),
    	       pyfits.Column(name='S_QUAN',
    			     format='1E',
    			     array=s_quan),
    	       pyfits.Column(name='LEV_DEG',
    			     format='1J',
    			     array=lev_deg),
    	       pyfits.Column(name='PHOT_TYPE',
    			     format='1J',
    			     array=phot_type),
    	       pyfits.Column(name='PHOT_PAR',
    			     format='20E',
    			     array=phot_par),
    	       pyfits.Column(name='ENERGY_REF',
    			     format='20A',
    			     array=energy_ref),
    	       pyfits.Column(name='PHOT_REF',
    			     format='20A',
    			     array=phot_ref)]
    	       ))

    hdu1.header.update('XTENSION', hdu1.header['XTENSION'],
    		    comment='Written by '+os.environ['USER']+now.strftime('%a %Y-%m-%d %H:%M:%S')+ 'UTC')
    hdu1.header.update('EXTNAME', ionstr, comment='Ion Name', before="TTYPE1")
    hdu1.header.update('HDUCLASS', 'ATOMIC',
    		    comment='Atomic Data', before="TTYPE1")
    hdu1.header.update('HDUCLAS1', 'E_LEVEL',
    		    comment='Energy level tables', before="TTYPE1")
    hdu1.header.update('ELEMENT', z,
    		    comment='Numer of protons in element', before="TTYPE1")
    hdu1.header.update('ION_STAT', q,
    		    comment='ion state (0 = neutral)', before="TTYPE1")
    hdu1.header.update('ION_NAME', ionstr, comment='Ion Name', before="TTYPE1")
    hdu1.header.update('N_LEVELS',len(energy) ,
    		     comment='Number of energy levels', before="TTYPE1")
    hdu1.header.update('HDUVERS1', '1.1.0',
    		     comment='Version of datafile', before="TTYPE1")
    if ionpot is not None:
      hdu1.header.update('ION_POT',ionpot,
             comment='Ionization potential (eV)',before="TTYPE1")

    for icmt in comments:
      hdu1.header.add_comment(icmt)

    # combine hdus
    hdulist = pyfits.HDUList([hdu0,hdu1])

    # write out file (overwrite any existing file)
    if os.path.isfile(filename): os.remove(filename)
    hdulist.writeto(filename, checksum=True)

# ------------------------------------------------------------------------------

  def write_xstar(self,filename,subset=None,precision=5):
    """ 
    Write xstar type 6 file. Must have ionization potential and nuclear charge
    specified.
    """
    if self.ionization_threshold() == None:
      raise ErrorRAP('must define ionization threshold')
    if self.nuclear_charge() == None:
      raise ErrorRAP('must define nuclear charge')
    indices=self.indices()
    if subset is not None: indices=subset

    zeff=self.nuclear_charge()
    ne=self.num_electrons()
    ionno=xstar_ionencode(zeff,ne)

    aform='%'+str(8+precision)+'.'+str(precision)+'E'
    iform='%6i'
    prefix=iform%6+iform%13+iform%0+iform%4+iform%6
    fout=open(filename,'w')
    for i in indices:
      conf=self.configuration(i)
      e=self.energy(i)*RYDBERG
      wgt=float(self.weight(i))
      ip=self.ionization_threshold()*RYDBERG
      np=-1.0
      if ip > 0.0: 
        z1=float(zeff-ne+1)
        de=abs(self.ionization_threshold()-self.energy(i))
        np=z1/sqrt(de)
      nval=0
      if conf is not None: nval=rlf_ssdecode(conf.valence_subshell())[0]
      (g,l)=rlf_termdecode(self.term(i))
      if g is None: 
        g,l=0,-1
      else:
        par=conf.parity()
        if par == 1: g=-g
      label=self.label(i,spacer='.')
      llen=len(label)
      line=prefix+iform%llen+' '+aform%e+aform%wgt+aform%np+aform%ip+' '+\
                  iform%nval+iform%g+iform%l+iform%zeff+iform%i+\
                  iform%ionno+' '+label+' %'
      fout.write(line+'\n')
    fout.close()

################################################################################

class rap_jconfig:
  """ Stores a single configuration in JJ-coupling. """
  def __init__(self,conforne,ne=None,ssocc=None):
    if isinstance(conforne,int):
      self.v_conf=self._baseconfig(conforne)
    else:
      self.v_conf=self._format_config(conforne,ne=ne,ssocc=ssocc)

# ------------------------------------------------------------------------------

  def __repr__(self):
    return self.label()

# ------------------------------------------------------------------------------

  def __eq__(self,conf):
    tconf=conf
    if not isinstance(tconf,rap_jconfig): tconf=rap_jconfig(tconf)
    if repr(self) == repr(tconf): return True
    return False

# ------------------------------------------------------------------------------

  def __ne__(self,conf):
    tconf=conf
    if not isinstance(tconf,rap_jconfig): tconf=rap_jconfig(tconf)
    if repr(self) != repr(tconf): return True
    return False

# ------------------------------------------------------------------------------

  def __cmp__(self,conf):
    return rlf_jconfsorter(self,conf)

# ------------------------------------------------------------------------------

  def __hash__(self):
    return hash(self.label())

# ------------------------------------------------------------------------------

  def __len__(self):
    """ Return length of configuration label. """
    return len(self.label())

# ------------------------------------------------------------------------------

  def __add__(self,other):
    return self.label()+other

# ------------------------------------------------------------------------------

  def __radd__(self,other):
    return other+self.label()

# ------------------------------------------------------------------------------

  def __getitem__(self,idx):
    return self.occupancy(idx)

# ------------------------------------------------------------------------------

  def __iter__(self):
    return self.subshells().__iter__()

# ------------------------------------------------------------------------------

  def _baseconfig(self,ne):
    """ Fill sub-shells in order until number of electrons is reached. """
    netot=ne
    out={}
    i=1
    while netot > 0:
      w=rlf_jssoccupancy(i)
      out[i]=min(w,netot)
      netot-=w
      i+=1
    return out

# ------------------------------------------------------------------------------

  def _format_config(self,inconf,ne=None,ssocc=None):
    """ Attempt to understand input configuration. """

    # configuration must be a string
    if isinstance(inconf,rap_jconfig):
      return inconf
    if not isinstance(inconf,str): 
      msg='configuration must be a string; conf: %s'%inconf
      raise ErrorRAP(msg)

    conf=inconf.lower()
    out={}

    # standard format
    tmp=re.findall(r'(\d+[a-z][+-])(\d+)',conf)
    if len(tmp) == 0:
      msg='cannot understand configuration string'
      raise ErrorRAP(msg)
    for (ss,tw) in tmp:
      w=int(tw)
      wmax=rlf_jssoccupancy(ss)
      if w > wmax:
        msg='occupancy exceeded for %s'%ss
        raise ErrorRAP(msg)
      iss=rlf_jssindex(ss)
      out[iss]=w

    # use subshell occupancies to figure out unspecified occupancies
    if ssocc is not None:
      for ss,ssw in ssocc.iteritems():
        (n,l)=rlf_ssdecode(ss)
        jss2=rlf_jssencode(n,l,+1,index=True)
        if l == 0:
          if jss2 not in out: out[jss2]=ssocc[ss]
          if out[jss2] != ssw:
            msg='configuration string incompatible with given sub-shell '+\
                'occupancy (%s)'%ss
            raise ErrorRAP(msg)
          continue
        jss1=rlf_jssencode(n,l,-1,index=True)

        if jss1 in out: ssw-=out[jss1]
        if jss2 in out: ssw-=out[jss2]
        if jss2 not in out:     # this subshell holds more electrons than jss1
          out[jss2]=0
          if ssw >= rlf_jssoccupancy(jss2):
            out[jss2]=rlf_jssoccupancy(jss2)
            ssw-=out[jss2]
        if jss1 not in out:
          out[jss1]=0
          if ssw > 0:
            out[jss1]=rlf_jssoccupancy(jss1)
            ssw-=out[jss1]
        if ssw != 0:
          msg='given sub-shell occupacy incompatible with given '+\
              'configuration: %s -- %s'%(inconf,ss)
          raise ErrorRAP(msg)          

    # fill in closed sub-shells
    issmin=min(out.keys())
    for iss in range(1,issmin): out[iss]=rlf_jssoccupancy(iss)

    # complete electron count
    if ne is not None:
      netot=sum(out.values())
      if netot > ne:
        msg='configuration has too many electrons, %s'%inconf
        raise ErrorRAP(msg)
      elif netot < ne:
        neleft=ne-netot
        ssmax=max(out.keys())
        for iss in range(1,ssmax):
          if neleft <= 0: break
          if iss in out: continue
          w=rlf_jssoccupancy(iss)
          if w <= neleft:
            out[iss]=w
            neleft-=w
            
        nechk=sum(out.values())
        if nechk != ne:
          msg='could not match given number of electrons'
          raise ErrorRAP(msg)

    # remove empty final subshells
    while 1:
      sss=sorted(out.keys())
      if out[sss[-1]] > 0: break
      del out[sss[-1]]

    return out

# ------------------------------------------------------------------------------

  def copy(self):
    """ Return a copy. """
    return copy.deepcopy(self)

# ------------------------------------------------------------------------------

  def label(self,showclosed=False):
    if self.num_electrons() == 0: return ''
    out=[]
    if showclosed:
      for ss in self.subshells(): out.append(ss+str(self.occupancy(ss)))
    else:
      imax=max(self.v_conf.keys())
      icore=1
      while icore in self.v_conf and \
            self.v_conf[icore] == rlf_jssoccupancy(icore):
        icore+=1
      if icore > imax: 
        icore=imax
      elif icore not in self.v_conf:
        out.append(rlf_jsubshell(icore)+'0')
      for i in self.subshells(indices=True):
        w=self.occupancy(i)
        if i < icore: continue
        out.append(rlf_jsubshell(i)+str(w))
    out=' '.join(out)
    return out

# ------------------------------------------------------------------------------

  def num_electrons(self):
    """ Return total number of electrons. """
    return sum(self.v_conf.values())

# ------------------------------------------------------------------------------

  def subshells(self,indices=False):
    """ Return sorted list of occupied sub-shells. """
    ss=self.v_conf.keys()
    ss.sort()
    if indices: return ss
    return map(rlf_jsubshell,ss)

# ------------------------------------------------------------------------------

  def occupancy(self,ss):
    """ 
    Return occupancy of sub-shell, ss (which can be index or string form. 
    """
    sss=ss
    if not isinstance(sss,int): sss=rlf_jssindex(sss)
    if sss not in self.v_conf: return 0
    return self.v_conf[sss]

# ------------------------------------------------------------------------------

  def parity(self):
    """ Return parity of configuration. """
    p=1
    for i in self.subshells():
      (n,l,s)=rlf_jssdecode(i)
      if l%2 == 0: continue
      p*=(-1)**self.occupancy(i)
    return (1-p)/2

# ------------------------------------------------------------------------------

  def weight(self):
    """ Return statistical weight of configuration. """
    wgt=1
    for ss in self.subshells():
      w=self.occupancy(ss)
      wmax=rlf_jssoccupancy(ss)
      wgt*=rlf_binom(wmax,w)
    return int(wgt)

# ------------------------------------------------------------------------------

  def valence_subshell(self):
    """ Return valence sub-shell. """
    tss=max(self.subshells(indices=True))
    return rlf_jsubshell(tss)

# ------------------------------------------------------------------------------

  def pop_electron(self,ss=None):
    """ 
    Return a copy of the configuration where an electron from the given 
    sub-shell has been removed (last sub-shell by default).
    """
    out=copy.deepcopy(self)
    tss=ss
    if tss is None: 
      tss=max(out.subshells(indices=True))
    else:
      if not isinstance(tss,int): tss=rlf_jssindex(tss)
    if tss not in out.subshells(indices=True): 
      msg='specified sub-shell not present: %s'%ss
      raise ErrorRAP(msg)
    out.v_conf[tss]-=1
    if out.v_conf[tss] == 0: del out.v_conf[tss]
    return out

# ------------------------------------------------------------------------------

  def add_electron(self,ss=None):
    """ 
    Return a copy of the configuration where an electron has been added to the
    given sub-shell (valence sub-shell by default).
    """
    out=copy.deepcopy(self)
    tss=ss
    if tss is None: 
      tss=max(out.subshells(indices=True))
    else:
      if not isinstance(tss,int): tss=rlf_jssindex(tss)
    if tss not in out.subshells(indices=True):
      out.v_conf[tss]=1
    else:
      if out.v_conf[tss] < rlf_jssoccupancy(tss):
        out.v_conf[tss]+=1
      else:
        if ss is None:
          out.v_conf[tss+1]=1
        else:
          msg='specified sub-shell is full: %s'%ss
          raise ErrorRAP(msg)
    return out

# ------------------------------------------------------------------------------

  def difference(self,conf):
    """ 
    Return dictionary of sub-shells and difference of occupation numbers.
    The sign of the difference is according to self - conf.
    """
    sss=self.subshells(indices=True)
    osss=conf.subshells(indices=True)
    out={}
    for ss in sss:
      w=self.occupancy(ss)
      if ss not in osss:
        out[rlf_jsubshell(ss)]=w
      else:
        dw=w-conf.occupancy(ss)
        if dw != 0: out[rlf_jsubshell(ss)]=dw
    for ss in osss:
      if ss in sss: continue
      out[rlf_jsubshell(ss)]=-conf.occupancy(ss)
    return out

# ------------------------------------------------------------------------------

  def to_config(self):
    """ 
    Convert configuration to rap_config form.
    """
    tmp={}
    for jss in self:
      ss=jss[:-1]
      if ss not in tmp: tmp[ss]=0
      tmp[ss]+=self[jss]
    keys=tmp.keys()
    keys.sort(rlf_sssorter)
    conf=[ss+str(tmp[ss]) for ss in keys]
    conf=' '.join(conf)
    return rap_config(conf)

